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Abstract 

We present a thorough study of the static properties of 2D models of spin-ice 
type on the square lattice or, in other words, the sixteen-vertex model. We use 
extensive Monte Carlo simulations to determine the phase diagram and critical 
properties of the finite dimensional system. We put forward a suitable mean-field 
approximation, by defining the model on carefully chosen trees. We employ the 
cavity (Bethe-Peierls) method to derive self-consistent equations, the fixed points 
of which yield the equilibrium properties of the model on the tree-like graph. We 
compare mean- field and finite dimensional results. We discuss our findings in the 
context of experiments in artificial two dimensional spin ice. 
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1 Introduction 



Many interesting classes of classical and quantum magnetic systems are highly con- 
strained. In particular, geometric constraints lead to frustration and the impossibility 
of satisfying all competing interactions simultaneously. In most cases this phenomenon 
gives rise to the existence of highly degenerate ground states [HE], often associated with 
an excess entropy at zero temperature. Anti-ferromagnets on a triangular lattice and 
spin-ice samples are materials of this kind. 

In conventional spin-ice [3J (see [U [5] for reviews) magnetic ions form a tetrahedral 
structure in 3D, i.e. a pyrochlore lattice. This is the case, for instance, of the Dy +3 ions 
in the Dy2Ti207 compound. Since the /-electron spins have a large magnetic moment, 
they can be taken as classical variables at low enough temperature, say, T < 10 K, behav- 
ing as Ising doublets, forced to point along the axes joining the centers of the tetrahedra 
shared by the considered spin. The origin of geometric frustration in these systems is 
twofold: it arises from the non-colinearity of the crystal field and the effective ferromag- 
netic exchange, and from long-range couplings between the spins [5]. Since within this 
Ising formulation the long ranged dipolar interactions are almost perfectly screened at low 
temperature [61 [7] , in a simplified description only short-range ferromagnetic exchanges 
are retained [3J. Thus, topological frustration arises from the fact that the Ising axes in 
the unit cell are fixed and different, forced to point towards the centers of neighboring 
tetrahedra. The configurations that minimize the energy of each tetrahedron are the six 
states with two-in and two-out pointing spins. 

The ground state of the system is more easily visualized by realizing that each tetra- 
hedron in 3D can be considered as a vertex taking one out of six possible configurations 
on a lattice of coordination four. With this mapping the magnetic problem described 
above becomes equivalent to the model of water ice introduced in [SJ In this context, 
the entropy of the ground state satisfying all the local (two-in - two-out) ice constraints, 
with all six vertices taken with the same statistical weight, was estimated by Pauling with 
a simple counting argument [9]. Pauling's result is very close to the earlier measurements 
performed by Giauque and Stout [TU] on water ice and to the ground-state entropy of the 
magnetic spin-ice sample measured in the late 90s [TT]. Experimentally, the Boltzmann 
weights of the six vertices can be tuned by applying pressure or magnetic fields along 
different crystallographic axes. Indeed, the extensions of Pauling's ice model to describe 
more general ferroelectric systems lead to 'ice-type models' [12]. We will come back to 
this issue below. 

The local constraint leads to many peculiar features that have been studied experi- 
mentally, numerically, and analytically. In the ground state, the total spin surrounding 
a lattice point is conserved and constrained to vanish according to the two-in - two-out 
rule. This fact has been interpreted as a zero- divergence condition on an emergent vector 
field Q2]. Spins are regarded as fluxes and, quite naturally, an effective fluctuating 
electromagnetism emerges, where each equilibrium configuration is made of closed loops 
of flux. This analogy can be used to derive power-law decaying spatial correlations of 
the spins [T4"| [15] . with a parameter dependent exponent, that were recently observed ex- 
perimentally via neutron scattering measurements [16] . Such criticality of the disordered 
ground state, also called spin-liquid or Coulomb phase, had been observed numerically 
first [T7] (see [TS] for a more general discussion). A detailed description of the Coulomb 
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phase can be found in [19]. 

Thermal (or other) fluctuations are expected to generate defects, in the form of ice- 
breaking rule vertices. In the electromagnetic analogy described above a defect cor- 
responds to a charge, defined as the number of outgoing arrows minus the number of 
ingoing arrows. As such, a tetrahedra with three-in and one-out spins contains a neg- 
ative charge of magnitude two and the reversed configuration a positive charge also of 
magnitude two. In turn, the four-in cells carry a charge minus four and the four-out 
ones a charge plus four. Such vertices should be present in the samples under adequate 
conditions. The possibility of observing magnetic monopoles and Dirac strings as be- 
ing associated to defects has been proposed by Castelnovo et al. [2U] and investigated 
experimentally by a number of groups JTHJ I2T] . 

Spin ice can be projected onto 2D Kagome planes by applying specially chosen mag- 
netic fields or pressure. Recently, the interest in spin-ice physics has been boosted by 
the advent of artificial samples [22] on simpler square lattices. These are stable at room 
temperature and have magnetic moments that are large enough to be easily observed 
in the lab, thus giving access to the micro-states that can be directly visualized with 
microscopy. Following the same line of reasoning exposed in the previous paragraphs, 
such 2D ice-type systems should be modeled by the complete sixteen-vertex model on a 
square lattice, where all kind of vertices are allowed. 

The exact solution of the 2D ice model (i.e., the six- vertex model) [23], and some 
generalizations of it in which a different statistical weight is given to the six allowed 
vertices [2H |25l [26] , was obtained by Lieb and Sutherland in the late 60s using the 
transfer matrix technique with the Bethe Ansatz. Soon after, Baxter developed a more 
powerful method to treat the generic six- and eight- vertex models [27] and founded in this 
way the theory of integrable systems (in the eight- vertex model vertices with four in-going 
or four out-going arrows are allowed). The equilibrium phase diagrams of the six- and 
eight- vertex models are very rich: depending on the Boltzmann weight of the vertices the 
systems can be found in a variety of different phases, such as a quasi long-range ordered 
spin liquid phase (SL), two ferromagnetic (FM) and one (or two) antiferromagnetic (AF) 
phases, separated by different types of transition surfaces in the phase diagram. In the 
six-vertex case the SL phase is critical similarly to what is observed in the 3D Coulomb 
phase. 

From a theoretical perspective integrable vertex models are of notable interest. Their 
static properties can be mapped into spin models with many-body interactions, loop 
models, three-coloring problems, random tilings, surface growth, alternated sign matrices 
and quantum spin chains. A comprehensive discussion of these mappings goes beyond 
the scope of the present work. For a review the interested reader may consult Ref . [25] . 
Some comments on these alternative representations, when useful, will be made in the 
text. 

Much less is known about the equilibrium (and dynamic [291 EH]) properties of the 
sixteen vertex model in two and three dimensions. As the experimental interest in classical 
frustrated magnets of spin-ice type is now focused on understanding the nature and the 
role of defects and their effects on the samples' macroscopic properties, it is worth to 
complete the analysis of such a generic model. The special experimental simplicity of two 
dimensional samples suggests to start from the 2D case. Moreover, it is worth trying to 
extend at least part of the very powerful analytic machinery developed for the integrable 
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six- and eight- vertex models to the (non-integrable) complete sixteen-vertex models, 
where all defects are allowed. 

In this paper we describe the equilibrium properties of spin-ice systems with defects. 
We proceed in two directions. On the one hand, we study the static properties of the 
sixteen-vertex model on a square lattice by using Monte Carlo simulations in equilibrium. 
For simplicity, and in order to keep the model close to experimental realizations, we 
assume that the statistical weight of the defects remains relatively small (as will be 
explained carefully in the text). We establish its phase diagram and its critical properties, 
that we compare to the ones of the integrable models. On the other hand, we introduce 
a suitable mean-field approximation, defined on a well-chosen tree of vertices, and we 
adapt the cavity (Bethe-Peierls) method to derive self-consistent equations on such trees, 
the fixed points of which yield the exact solution of the mean-field model. This method 
allows us to describe all expected phases and to unveil some of their properties, such 
as the presence of anisotropic equilibrium fluctuations in the symmetry broken phases. 
We discuss the range of validity of the mean-field approximation, and we compare the 
analytical solution to the numerical results for the finite dimensional system. In the 
conclusions we summarize our results and we briefly discuss some experimental features 
that we can describe with our methods |31j . 

2 Vertex models 

In vertex models the degrees of freedom (Ising spins, q— valued Potts variables, etc.) sit 
on the edges of a lattice. The interactions take place on the vertices and involve the spins 
of the neighboring edges. In this Section we recall the definition and main equilibrium 
properties of Ising-like vertex models in two dimensions defined on a square lattice. 

2.1 Definitions 

We focus on an L x L square lattice with periodic boundary conditions. We label the 
coordinates of the lattice sites by (m,n). This lattice is bipartite, namely, it can be 
partitioned in two sub-lattices A± and A 2 such that the sites having m + n even belong 
to Ai, those having m + n odd belong to A2, and each edge connects a site in A\ to one 
in Ai- The degrees of freedom sit on the links between two sites or, in other words, on 
the "medial" lattice whose sites are placed on the midpoints of the links of the original 
lattice. The midpoints are hence labeled by (m + 1/2, n) and (m, n + 1/2). Thus, in the 
models we consider the degrees of freedom are arrows aligned along the edges of a square 
lattice, which can be naturally mapped into Ising spins, say s m +i/2,n = ±1- Without lose 
of generality, we choose a convention such that s — +1 corresponds to an arrow pointing 
in the right or up direction, depending on the orientation of the link, and conversely 
s = — 1 corresponds to arrows pointing down or left. 

2.2 The six- vertex model 

In the six- vertex model (i.e., 2D spin-ice) arrows (or Ising spins) sit on the edges of a 
coordination four square lattice and they are constrained to satisfy the two-in two-out 
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rule [H [9] . Each node on the lattice has four spins attached to it with two possible direc- 
tions, as shown in Fig. [TJ A charge can be attributed to each single vertex configuration, 
simply defined as the number of out-going minus the number of in-coming arrows. Ac- 
cordingly, the six-vertex model vertices have zero charge. (Note that the charge is not 
the sum of the spins attached to a vertex, such a total spin will be defined and used in 
Sec. H) 

Although in the initial modeling of spin-ice all such vertices were equivalent, the 
model was later generalized to describe ferroelectric systems by giving different statistical 
weights to different vertices: oc e~ l3ek with the energy of each of the k = 1, . . . , 6 
vertices. /3 = l/{ksT) is the inverse temperature. Spin reversal symmetry naturally 
imposes that W\ = W2 = a for the first pair of ferromagnetic (FM) vertices, = = b 
for the second pair of FM vertices, and w§ = w§ = c for the anti-ferromagnetic (AF) 
ones, see Fig. [U We have here introduced the conventional names a, b, and c of the three 
fugacities corresponding to the Boltzmann weight of the three kinds of vertices. In the 
literature it is customary to parametrize the phase diagram and equilibrium properties 
in terms of a/c and b/c. This is the choice we also make in this paper. Particular cases 
include: the F model of anti-ferroelectrics, in which the energy of the AF vertices is set 
to zero and all other ones are taken to be equal and positive, i.e. c> a = b [32]; the KDP 
model of ferroelectrics, in which the energy of a pair of FM vertices is set to zero and all 
other ones are equal and positive, e.g. a > b = c [T2]; the Ice model, in which the energies 
of all vertices are equal, i.e. a = b = c [23J. It is important to note, however, that in the 
context of experiments in artificial spin-ice type samples, vertex energies are fixed and the 
control parameter is the temperature. In [31] we used this alternative parametrization 
and we compared the model predictions to experimental observations [331 
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Figure 1: (Color online.) The zero-charge six- vertices of the six- vertex model. v\, V2, V3 and 
are FM vertices while V5 and vq are AF ones. Vertices are grouped in spin- reversal symmetric 
pairs. 

The six-vertex model is integrable. The free-energy density of the model with a = 
b = c and periodic boundary conditions was computed by Lieb in the late 60s with the 
transfer matrix technique and the Bethe Ansatz to solve the eigenvalue problem [25] . 
The method was then extended to calculate the free-energy density of the F model [21], 
the KDP model [25], and also models with generic values of a, b, and c, and periodic 
boundary conditions [26] , and the same general case with antisymmetric [35] and domain 
wall boundary conditions [36] • The effect of the boundary conditions is indeed very subtle 
in these systems as some thermodynamic properties depend upon them [37] contrary to 
what happens in usual short-range statistical physics models. Very powerful analytic 
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methods such as the Yang-Baxter equations have been employed in the analysis of these 
integrable systems [27]. 

An order parameter that allows one to characterize the different phases is the total 
direct and staggered magnetization per spin 

(Mt) = \{{\ml\) + <K|» (1) 
with the horizontal and vertical fluctuating components given by 

L 2 m X ± = S m+1 / 2 ,n =t •S m+ i/2,n ) (2) 

(m,n)sAi (m,n)sA2 

L 2 m V ± = S m ,n+l/2 i S m>rl _|_i/2 • (3) 

(m,n)gAi (m,n)6j42 

The angular brackets (...) denote here, and throughtout the rest of this article, the 
statistical average. A finite value of (M + ) indicates the spontaneous breaking of Z 2 
symmetry, while a finite value of (M_) corresponds to the spontaneous breaking of Z 2 
and translational symmetry. 




Figure 2: (Color online.) The phase diagram of the six- (red solid lines) and eight- (dashed 
lines) vertex models. For the eight- vertex model the curves correspond to the projection on the 
d = plane. Only for d = (six-vertex model) the PM region becomes an SL (Coulomb) phase. 



The four equilibrium phases are classified by the anisotropy parameter 

a 2 + b 2 - c 2 



(4) 



2ab 

and they are the following. 

a- Ferromagnetic (a-FM) phase: A 6 > 1; i.e. a > b + c. Vertices v± and v 2 are favored. 
Spin reversal symmetry is broken. The lowest energy state in the full FM phase is doubly 
degenerate: either all arrows point up and right or down and left [i.e. M + = 1, with M + 
the magnetization density defined in eq. ([!])]. In this phase the system is frozen as the 
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only possible excitations involve a number of degrees of freedom of the order of L. In all 
this phase the exact free-energy per vertex is given by [27] 



/fm — . (5) 

At a = b + c (A 6 = 1) the system experiences a frozen-to-critical phase transition from 
the frozen FM to a disordered (D) or spin liquid (SL) phase that we discuss below. 

b- Ferromagnetic (b-FM) phase: A$ > 1; i.e. b > a + c. This phase is equivalent to the 
previous one by replacing a- by 6-vertices. The free-energy is /fm = £3 and the phase 
transition towards the SL phase is also of frozen-to-critical type. 

Spin liquid (SL) phase: — 1 < Aq < 1; i.e. a < b + c, b < a + c and c < a + b. In this 
phase the averaged magnetization is zero, (M±) = 0, and one could expect the system to 
be a conventional paramagnet (PM). However, the ice constraints are strong enough to 
prevent the full decorrelation of the spins even at finite temperature. The system is in a 
quasi long-range ordered phase with an infinite correlation length. At c = a + b there is a 
Kosterlitz-Thouless (KT) phase transition between this critical phase and an AF phase 
with staggered order that is discussed below. Some particular points in parameter space 
belong to the SL phase as the spin-ice point a = b = c for which Aq = 1/2. At this 
special point the ground state is macroscopically degenerate giving rise to the residual 
entropy [23J 

S/N = 3/2 In 4/3 (6) 

with N = L 2 the number of vertices in the sample. 

The exact solution found by Baxter yields the free-energy density as a function of 
the parameters through a number of integral equations [27]. In Sec. H] we evaluate it 
numerically and we compare it to the outcome of the Bethe approximation. Close to the 
FM transitions the free-energy density can be approximated by 

/ SL ^ max(ei, e 3 ) - ^k B T {^~^- ~ A = max(ei, e 3 ) - ^k B T t 2 ~ a , (7) 

with t being the reduced distance from criticality, t = (6 + c)/a — 1, and a an exponent 
that plays the role of the one of the heat-capacity and takes the value a = 1 here. The 
first derivative of /sl with respect to the distance from the transition t shows a step 
discontinuity at the SL-FM transition as it would in a first-order phase transition, even 
though the FM phase is frozen. This corresponds to a critical-to-frozen phase transition. 

Antiferromagnetic (AF) phase: Aq < —1; i.e. c > a + b. Vertices v$ and v§ are favored. 
The ground state is doubly degenerate, corresponding to the configurations M_ = ±1. 
The staggered order is not frozen, due to thermal fluctuations. This is confirmed by the 
exact expression of the staggered magnetization found by Baxter [3BJI3H]. The free-energy 
has an essential singularity at the critical temperature (towards the SL phase) 

/af ^ e -cstM t (8) 

with est being a constant and t = (a + 6)/c — 1 the distance from criticality, as it is typical 
for an infinite order phase transition. 

The transition lines are straight lines (given by Aq = 1 for the SL-FM and Aq = — 1 
for the SL-AF) and they are shown in Fig. |2] as solid (red) lines. The dashed line along 
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the diagonal rapresents the range of variation of the parameters in the F model. The 
horizontal dashed line is the one of the KDP model. The intersection of this two lines 
corresponds to the ice-model. Although the transitions are not of second order, critical 
exponents have been defined and are given in the first column of Table [TJ for the SL-FM 
transition. The exponent a is taken from the expansion of the free-energy close to the 
transition, see eqs. (0) and (jHJ). The ratios 7 = 7/1/, j3 — j3/u and = (2 — a)jv are 
defined using the (divergent in the SL phase) correlation length £ instead of t as the 
scaling variable [10] • The values of 7, $ and <fi are the same for the two transitions and, 
interestingly enough, coincide with the ones of the 2D Ising model recalled in the last 
column of these tables. They also satisfy the usual scaling relations. The exponents j3, 7 
and v are extrapolated from the eight- vertex model ones as explained in the next section. 

2.3 The eight-vertex model 

The eight-vertex model is a generalization of the six-vertex model, first introduced to 
get rid of some of its very unconventional properties due to the hard ice-rule constraint 
(frozen FM state, quasi long-range order at infinite temperature, etc.) [H] 112] . In this 
model the allowed local configurations are those for which each vertex is surrounded by 
an even number of arrows pointing in or out, resulting in the addittion of the two vertices 
with four ingoing and four outgoing arrows shown in Fig. [3] to the ones in Fig. [TJ The 
eight-vertex model can still be mapped into a loop model on the lattice in such a way 
that the integrability property is preserved. It was first solved by Baxter in the zero-field 
case (i.e., with Z2 symmetry) [381 ES]. In order to do that he introduced the celebrated 
Star- Triangle relations (now called Yang-Baxter equations). 

Vi V 2 V 3 V4 V b Vq V 7 V 8 




a—coi—u>2 b—ojs—uj^ C—CO5—0JQ d—cuj—ujs 



Figure 3: (Color online.) The vertices in the eight-vertex model. The four-out and four-in 
vertices, vj and vg, with charge +4 and —4, respectively, have weight d. 

The phase diagram is characterized by the anisotropy parameter 

A 8 = a2 + fe2 ~ c2 -' 2 (9) 
8 2{ab + cd) 1 ' 

which becomes the six- vertex one when d = 0. This model sets into the following five 
phases depending on the weight of the vertices: 

a- Ferromagnetic phase (o-FM): A§ > 1 (a > b + c + d). Spin reversal symmetry is broken. 
This ordered phase is no longer frozen and (M+) < 1. 
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b- Ferromagnetic phase (6-FM): A§ > 1 (6 > a + c + d). This phase is equivalent to the 
previous one replacing a- by b- vertices. 

Paramagnetic phase (PM): —1 < A§ < 1 [a, b, c, d < (a + b + c + d)/2]. As soon as 
d > this phase is truly disordered, with a finite correlation length. The averaged 
magnetization vanishes (M±) = 0. 

c-Antiferromagnetic phase (c-AF): A§ < — 1 (c > a + b + d). Translational symmetry is 
broken. The configurations are dominated by c-vertices with an alternating pattern of 
vertices of type v 5 and v 6 with defects; (M_) < 1. 

d-Antiferromagnetic phase (d-AF): Ag < — 1 (d > a + b + c). Translational symmetry is 
broken. The configurations are dominated by ci-vertices, with an alternating pattern of 
vertices v-j and v 8 with defects. (M_) is also different from zero in this phase. This order 
parameter does not allow one to distinguish the d-AF from the c-AF. 

The transition lines are given by A 8 = 1 for the PM-FM ones and A 8 = — 1 for the 
PM-AF ones. The projection of the critical surfaces on the d = plane yields straight 
lines translated by d/c with respect to the ones of the six- vertex model, in the direction 
of enlarging the PM phase, as shown by the dashed lines in Fig. [2j 

The effect of the rf-vertices on the order of the different phase transitions is very 
strong. As soon as d > 0, the KT line between the c-AF and the SL phases becomes 
'stronger', i.e. second order, except at the intersection with the a = or b = planes 
when the transition is first order. On the contrary, the frozen-to-critical lines between 
the FMs and SL phases get "softer", i.e. second order, and become KT transitions on 
the a = and 6 = planes. Finally, the separation between the d-AF and disordered 
phases is second order for a, b, d > and it is of KT type on the a = and 6 = planes. 
As we will show in Sec. [31 this is consistent with our numerical results. 

The critical exponents can be found from the analysis of the free-energy density close 
to the transition planes and they depend explicitly on the weights of the vertices via the 
parameter tan(/x/2) = \Jcd/ab [38] . The critical exponents for the PM-FM transitions in 
this model are given in the second column of Table [TJ The values of /3, 7 and for the 
six-vertex model given in the first column are obtained as the limit of these exponents 
for d — > when ji — ir. 





six-vertex 


eight-vertex 
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1/16 


7T/(16/i) 


1/8 


7 
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7tt/(8m) 


7/4 
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1/2 
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Table 1: Exact critical exponents of the six- and eight- vertex model with tan(///2) = y 'cd/ab 
along the SL/PM-FM transition. In the limit d — > the parameter fi — > tt and the eight- vertex 
model exponents become the ones of the six-vertex model. The critical exponents of the 2D 
Ising model are recalled in the third column for comparison. 
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2.4 The sixteen vertex model 



The most general model obtained by removing the ice-rule is the sixteen- vertex model, in 
which no restriction is imposed on the value of the binary variables attached on each edge 
of the lattice, and all the 2 4 = 16 vertex configurations can occur. The (eight) three-in 
one-out and three-out one-in vertices that are added to the ones already discussed are 
shown in Fig. HI In order to preserve the Z2 symmetry (in absence of an external magnetic 
field that would break the rotational symmetry), the same statistical weight e is given 
to all these 'defects' with charge 2 and —2. In the figure, vertices are ordered in pairs of 
spin-reversed couples [v\q is the spin reversed of Vg and so on) and the difference with 
the following couples is a rotation by it (vn is equal to Vg apart from a 7r-rotation and so 
on). 



Vg 



V10 



Vn 



V12 



V13 



V14 Vi 5 



VlQ 



> o > 



< r Hh < r < r > r Hh 4 



A V A 

< > • > ><<-<> 



e=ojg=a;io = ...=wi6 



Figure 4: (Color online.) The eight three-in one-out (with charge —2) or three-out one-in (with 
charge +2) vertices that are included in the sixteen vertex model. We give them equal weight 
e. 

The new vertices naturally entail the existence of new phases. One can envisage the 
existence of a critical SL phase for a = b = c = d = ^10,12,14,16 = an d ^11^3^5 > as 
this new four- vertex model is equivalent to the dimer model solved by Kasteleyn [13]. It 
is quite easy to see that e-AF stripe order is also possible. For instance, one can build an 
ordered configuration with alternating lines of Vg and viq vertices, or another one with 
alternating columns of Vn and V12 vertices. Phases of this kind should appear if one favors 
one pair of spin-reversed related vertices by giving them a higher weight than the others, 
and the transition to this phase should be continuous, since local fluctuations made by 
elementary loops around a square plaquette are allowed. 

The sixteen- vertex model loses the integrability properties j28j HU H5] . Nonetheless, 
some exact results are available for a few special sets of parameters. The equivalent 
classical Ising model has only nearest and next-nearest neighbor two-body interactions 
when e 4 = abed fZE\ HI]. In the c-AF sector this condition leads to the generalised F 
model defined bjc=l,a = b<l,d = a u and e = a v , with the constraint Av — u + 2. 
The model has been solved for the special cases: (i) v — 1 and u = 2 (i.e. e = a and 
d = a 2 ), the associated spin model simplifies into an AF Ising model with only nearest 
neighbor interactions. This model is known to exhibit a second-order phase transition^ 
with a logarithmic divergence of the specific heat (a = 0). (ii) Using a different approach 

1 The transition occurs at e/fcsT c = 21n(\/2 + 1) ~ 0.567. Using our conjectured Ai 6 (denned below) 
we get a critical temperature £//cbT c w 0.607. 
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it has been showed that for v — > oo and u — 2 (i.e. e = and d = a 2 ) the system also 
exhibits a second-order phase transition in the same universality class as (i). Notice that 
the exactly solvable F model is recovered in the limit v — > oo and u — > oo. In the same 
way, in the a-FM sector this leads to the generalised KDP model [IB] by setting a — 1, 
b = c < 1, d = b u , e = b v and again 4v = u + 2. For v — 1 and u — 2 {i.e. e = b and 
<i = b 2 ) the system exhibits a second-order phase transition with the same properties of 
its c-AF analog discussed above. For v — > oo and u = 2 (i.e. e = and d = b 2 ) the 
system also exhibits a second-order phase transition in the same universality class as the 
previous case. 

Since the phase diagram of the generic sixteen-vertex model is rather complex, and 
our aim is to consider e and d vertices as defects having a relative small statistical weight, 
in this work we will focus on the effect of the presence of e vertices on the phases and the 
phase transition lines described in Fig. |2] 

2.5 Numerical simulations 

The numerical analysis of the equilibrium properties of 2D vertex models has been re- 
stricted, so far, to the study of the six and eight-vertex cases. Single spin-flip updates 
break the six- and eight- vertex model constraints and cannot be used to generate new 
allowed configurations. Instead, as each spin configuration can be viewed as a non- 
intersecting (six-vertex) or intersecting (eight-vertex) loop configuration, stochastic non- 
local updates of the loops have been used to sample phase space [HI 0SJ. By imposing 
the correct probabilities all along the construction of non-local moves, cluster algorithms 
can be designed j5HJ S3] • Non-trivial issues as the effect of boundary conditions have been 
explored in this way [HJ [50] . 

Loop-algorithms, as usually presented in the context of Quantum Monte Carlo meth- 
ods, exploit the world-line representation of the partition function of a given quantum 
lattice model [51]. It is well known that the 2D six- and eight- vertex models are equiv- 
alent to the Heisenberg XXZ and XYZ quantum spin-1/2 chains, respectively |41[ [52] . 
It is then not surprising to find the same kind of loop-algorithms in the vertex models 
literature. A configuration in terms of bosonic world lines of the quantum spin chain 
in imaginary time can be one-to-one mapped into a vertex configuration on the square 
lattice, so that the same loop algorithm samples equivalently the configurations of both 
models. 

The loop algorithms could be modified to include three-in - one-out and one-in - 
three-out defects for the study of spin-ice systems [531 [M]- However, the simultaneous 
inclusion of four-in and four-out defects makes this algorithm inefficient compared to a 
Monte Carlo algorithm with local updates. For this reason, we will use local moves in 
our numerical studies, as we explain in the next Section. 

3 Numerical study of the sixteen vertex model 

In this Section we summarize the results obtained for the sixteen-vertex model using 
Monte Carlo (MC) simulations. We first explain the numerical algorithm. Then, we 
discuss the phase diagram and the critical properties of the model. All our results are for 
a square lattice with linear size L and periodic boundary conditions. 
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3.1 Numerical methods 



We used two numerical methods to explore the equilibrium properties of the generic 
model; the Continuous time Monte Carlo (CTMC) method that that we briefly explain 
in Sec. 13.1.11 and in App. |A] and the Non-equilibrium relaxation method (NERM) that 
we equally briefly explain in Sec. 13.1.21 

3.1.1 Continuous time Monte Carlo 

The usual Metropolis algorithm, from now on called Fixed Step Monte Carlo algorithm 
(FSMC), is very inefficient to study the equilibrium properties of frustrated magnets. 
The dynamics freeze when d, e <C min(a,6, c) and the acceptance probability for most 
of the single-spin flip updates is extremely small. We implemented an algorithm which 
overcomes this difficulty, the Continuous Time Monte Carlo algorithm (CTMC) [55| 156]. 
This method is also known with different names: Bortz-Kalos-Lebowitz [55], n-fold way 
or kinetic Monte Carlo. The basic aim of this algorithm is to get rid of the time wasted 
due to a large number of rejections when the physics of the problem imposes a very 
small acceptance ratio. It is extremely useful for the study of the long time behavior of 
systems with complicated energy landscapes and a large number of metastable states. 
The main idea behind the method is to sample stochastically the time needed to update 
the system and then do it without rejections. In our case, we chose to use single spin 
updates such that the two vertices connected by the spin are updated in the move. Details 
on this algorithm are given in App. |A] Equilibrium can be achieved for relatively large 
samples. The finite size scaling analysis of the equilibrium Monte Carlo data yield the 
thermodynamic properties of the generic model. 

3.1.2 Non-equilibrium relaxation method 

The fact that dynamic scaling applies during relaxation at a critical point [57] suggested 
to use short-time dynamic measurements to extract equilibrium critical exponents with 
numerical methods [5E1 EE1 EH] • Magnetized, M° = M±(t = 0) ^ 0, and non- magnetized, 
M± = 0, configurations can be used as starting conditions and the critical relaxation 

M±{t) ~ t~ PI{vz) F{t Xo/z Ml) (10) 

where z is the dynamic critical exponent, and F(x) ~ x for i < 1 and F(x) —> est for 
x — > oo, can be used to extract either the critical parameters or the critical exponents. 
This expression is expected to hold for t l l z <C L and t x l z <C ^ eq with £ eq the equilibrium 
correlation length. 

3.2 Phase diagram and critical singularities 

In this subsection we present a selected set of results from our simulations, and we describe 
the kind of phases and critical properties found. The strategy to study the different phase 
transitions is the following. We first chose the relevant order parameter, (M + ) or (M_), 
to study FM or AF phases. From the finite size scaling analysis of the corresponding 
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fourth-order reduced Binder's cumulant 



K M± = 1 - - 9 K (tL^) (11) 

where t is the distance from the critical point, we extracted the critical exponent v. From 
the maximum of the magnetic susceptibility 

X± = L 2 [(M 2 ± ) - (M ± ) 2 ] ~ L^ x (tL^) (12) 

we extracted j/v, then 7 as v was already known. From the maximum of the specific 
heat 

C E = L~ 2 [(E 2 ) - (E) 2 ] ~ L a/ ^ c (tL 1/u ) (13) 

we extracted a/v, then a. Here, E = n^k with n\~ the number of vertices of type k 
and ek their energy. The direct measurement of (3 is difficult, we thus deduced it from 
the scaling relation (3 = |(2 — a — 7). Finally, we checked hyper- scaling, i.e. whether 
dv = 2 — a is satisfied by the exponent values obtained, that we summarize in Table |2] 
for the SL/PM-FM transition and two choices of parameters. 



3.2.1 The PM-FM transition 

In order to reduce the number of parameters in the problem we studied the PM-FM 
transition for the special choice d = e. 

As the direct magnetization density (M + ) defined in eqs. ([I])-© is the order parameter 
for the PM-FM transition in the six-vertex model, we study this quantity to investigate 
the fate of the FM phase in the sixteen-vertex model. In Fig. [5] we show (M+) as a 
function of a for b = 0.5 and three values of the fugacity d = e (all normalized by c). 
The data for the 2D model (shown with colored points) demonstrate that the presence of 
defects tends to disorder the system and, therefore, the extent of the PM phase is enlarged 
for increasing values of d — e. Moreover, the variation of the curves gets smoother for 
increasing values of d = e suggesting that the transitions are second order (instead of 
frozen-to-critical) in presence of defects. The data displayed with black points for the 
same parameters are the result of the mean-field analysis of the model, which will be 
discussed in Sec. HI In [2H] we suggested that the equilibrium phases of the sixteen vertex 
model with d = e could be characterized by a generalization of the anisotropy parameter 
of the eight- vertex model recalled in eq. (Q: 

a 2 + fe 2 -c 2 -(rf + 3e) 2 
16 2[ab + c(d + 3e)} ' 1 j 

In the same way as in the integrable cases, the proposal is that the PM phase corresponds 
to the region of the parameters' space where |Ai6| < 1, the FM phases corresponds to 
A i 6 > 1, and the AF ones to A 16 < —1. It follows that the projection of the FM-transition 
hyper-planes onto the (a/c, b/c) plane should then be parallel to the ones of the six- and 
eight- vertex models and given by a c = b + c + d + 3e (or equivalently b c = a + c + d + 3e). 
As shown in Fig. |5j the numerical results are well fitted by Eq. f!14p which is, however, 
not exact. As we will explain in detail in Sec. HJ our mean- field treatment of the model 
defined on a tree of single vertices predicts a similar shift of the transition lines given by 
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a c = b + c + d + 2e. The parameter capturing all transition lines is, within this analytic 
approach, 

ASV _ a 2 + b 2 -c 2 -d 2 + 2{a + b-c- d)e 
16 ~ 2(cd + ab + e(a + b + c + d + 2e)) ^ ^ 

(see Sec. H]for the technical details). For the more sophisticated tree made of 'plaquette' 
units the transition lines are not parallel to the ones of the six- and eight-vertex models 
and an analytic form of the anisotropy parameter, A^g, has not been found. Nevertheless, 
the transition lines can be computed numerically and their evaluation leads to the phase 
diagram depicted in Fig. [T9l 




Figure 5: (Color online.) Equilibrium magnetization density, (M+), of the sixteen vertex 
model for three different values of d = e (given in the key) and b = 0.5 as a function of a 
(vertices fugacities are normalized by c). The colored data points are the result of the numerical 
simulations of the 2D model for L = 40 while the black dots corresponds to the analytic solution 
of the model defined on the tree of plaquettes, as explained in Sec. HI 

Further evidence for the transition becoming second order comes from the analysis 
of the fourth-order cumulant defined in eq. ffTTj) . Raw data on such Binder's cumulant 
across the FM-PM transition were shown in [29] where we showed that they intersect at a 
single point, as expected in a second order phase transition. In Fig. Owe display raw data 
for d = e = 10~ 5 (a) and scaled data for d = e = 0.1 (b) as a function of t — (a — a c )/a c . 
In both cases b = 0.5 and, as above, we normalize all fugacities by c. From the analysis of 
the scaling properties we extract a c = 1.5 for d = e = 10 -5 and a c = 1.93 for d = e = 0.1. 
Sets of data for linear system sizes L = 10, 20, 30, 40, 50 are scaled quite satisfactorily 
by using \ jv — 1.65 ± 0.05 for the small d and 1/u = 1 zb 0.1 for the large value of d. 

In order to complete the analysis of this transition we studied the magnetic suscepti- 
bility p2|) associated to the direct magnetization M + and its finite size scaling. Figure [7] 
displays x+ f° r b = 0.5, d — e — 10~ 5 (a) and d — e — 0.1 (b), and five linear sizes, 
L = 10, 20, 30, 40, 50. The data collapse is very accurate and it allows us to extract the 
exponent y/V ~ 1.75 ± 0.02 in both cases. The study of the maximum of x+ displayed 
in the insets confirms this estimate for y/iA We repeated this analysis for other values of 
d — e and we found that in all cases critical scaling is rather well obeyed and, interestingly 
enough, 7/1/ is, within numerical accuracy, independent of d = e. This is similar to what 
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Figure 6: (Color online.) Analysis of the Binder fourth-order cumulant defined in eq. (jlip 
across the FM-PM transition in the sixteen vertex model, (a) Raw data for d = e = 10 -5 . (b) 
Scaling plot for d = e = 0.1. One extracts \/v = 1.65 ± 0.05 in case (a) and 1/v = 1 ± 0.1 in 
case (b) from this analysis. 



happens in the eight- vertex model where, as shown in Table (H this ratio is independent 
of the parameters. 

The ratio ajv is obtained from the finite size analysis of the specific heat (not shown) 
that is consistent with C max ~ L a l v (instead of L D for a first order phase transition). We 
found ajv ~ 1.30 ±0.06 for b = 0.5 and d = e = 10 -5 and a logarithmic divergence of the 
heat capacity, i.e. ajv »s for d = e = 0.1 (cf. Table [2]), although it is very difficult to 
distinguish numerically a logarithmic divergence from a power-law one with a very small 
exponent. 

The critical exponents extracted numerically at the second order phase transition with 
d = e > are compared to the ones of the six-vertex model and the 2D Ising model in 
Table [2j It is interesting to note that for very small value of d = e the exponents are rather 
close to the ones of the six-vertex model while for large value of d = e they approach 
the ones of the 2D Ising model. In terms of a Renormalization Group approach, this 
behavior suggests the existence of two fixed points, one in the d = e = plane describing 
the critical behavior of the six-vertex model, and another for d, e > 0, belonging to the 
2D Ising universality class. The first one might be unstable as soon as an infinitesimal 
amount of defects is allowed. A RG treatment of the model is required to confirm this 
guess. 

The critical exponents obtained from the numerical analysis depend on the fugacity 
d = e, namely, they vary along the transition lines, as it happens for the eight-vertex 
model. However, as shown in Table El the ratios of critical exponents 7, and <p do not 
depend on the choice of the parameters. As one could expect, the values we obtained 
for these ratios of critical exponents are equal to the ones of the eight-vertex model and 
the 2D Ising model, since these two latter models are special cases of the sixteen vertex 
model. 

The NERM yields values of the critical a that are in agreement (within numerical 
accuracy) with the ones found with the conventional analysis. We do not show this 
analysis here. 
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Figure 7: (Color online.) The magnetic susceptibility across the PM-FM transition for b = 0.5, 
d = e = (a) and d = e = 0.1 (b). Where t is the distance from the critical point measured 
as t = (a — a c )/a c with a c = 1.5 (a) and a c = 1.93 (b). From the finite size scaling of the 
maximum shown in the insets one extracts j/v ~ 1.75 ± 0.02 in both cases. 
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Table 2: Numerical values of the critical exponents at the FM-PM transition in the sixteen 
vertex model as compared to the ones in the six-vertex model (first column) and 2D Ising model 
(fourth column). The parameter fx, tan(///2) = W cd/ab, has been chosen to take the same value 
in the two MC columns, fi = ir. We did not include errorbars in the column corresponding to 
d = e = 0.1 as our determination of a is not precise enough to distinguish between a = (the 
value used to extract the remaining exponents) and a very small but non-vanishing value. 



3.2.2 The c-AF-PM transition 

We now focus on the transition between the c-AF and PM phases. For the six-vertex 
model this is a KT transition while for the eight- vertex model it is of second order as soon 
as d > 0. In this case we chose to work with d = 10~ 5 7^ e = 10~ 3 and with d = e = 10 -5 . 
We present data obtained with the NERM. 

Figure [H] shows the relaxation of the staggered averaged magnetization at b = c = 1 
and different values of a given in the caption. The power-law relaxation, typical of 
the critical point, is clearly identifiable from the figure. We extract the critical value 
a c = 0.46 ±0.01 for the c-AF-PM phase transition. Moreover, the data demonstrate that 
the PM phase is not of SL-type as soon as a finite density of defects is allowed. Indeed, 
the relaxation of M_ does not follow a power law in the PM phase, the decay being 
exponential for a > a c . Although this strategy gives a rather precise determination of 
a c , it is very hard to determine the value of the exponent p, and hence of z, with good 
precision as p is extremely sensitive to the choice of a c . 
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Figure 8: (Color online.) Non-equilibrium relaxation of the staggered magnetization from a 
fully ordered initial condition M®_ = 1 at different values of a = b, for c = 1, e = 10 -3 and 
d = 10 -5 . After a short transient the relaxation at the critical point follows a power law t~ p 
with p = (3/(vz). We identify such critical relaxation at a c = 0.46 ± 0.01. 

The standard analysis of the c-AF-PM transition is not as clean as for the FM-PM one. 
Figure [9] (a) shows the scaling plot of the Binder cumulant of the staggered magnetization 
for b = 0.5 and, in this case, d = e = 10 -5 . From it one extracts \ jv = 0.4 ± 0.05. The 
susceptibility fluctuates too much to draw any stringent conclusion about the exponent 
7. The analysis of the specific heat (not shown) suggests a logarithmic divergence a ~ 0. 
The dependence of the critical line on the fugacities is reasonably well described by Ai6- 

Both e and d vertices make the disordered phase be a conventional PM, and the transi- 
tions be second order. Although d does not break integrability while e does, their effect, 
in these respects, are similar. The fact that the critical exponents in the eight-vertex 
model depend on the fugacities is known from the exact solution. In the sixteen-vertex 
model, where integrability is lost, this is not the case. Unfortunately, our numerical anal- 
ysis does not allow us to draw stringent enough conclusions for the c-AF-PM transition, 
since it is very hard to get precise measurement of the observables close to the critical 
lines. Possibly, a renormalization group approach would be useful in this respect. 

4 Bethe-Peierls mean-field approximation: Vertex 
models on a tree 

In this Section we study the properties of the six-, eight- and sixteen-vertex models denned 
on the Bethe lattice by using the cavity method. First, we will consider a standard Bethe 
lattice of uniform connectivity C = 4 (Sec. I4.2.ip . In order to get more accurate results, 
we also study the model on a tree of plaquettes of 2 x 2 vertices (Sec. I4.2.2p . which 
account for local excitations and fluctuations induced on short scales by the presence of 
small loops. 
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Figure 9: (Color online) Study of the c-AF-PM transition for b = 0.5 and d = e = 10~ 5 . (a) 
The averaged staggered magnetization as a function of a for several system sizes given in the 
key. (b) Scaling plots across the c-AF-PM transition of the Binder cumulant of the staggered 
magnetization, t is the distance from the critical point t = (a — a c )/a c with a c = 0.5. The best 
scaling of data is obtained for 1/v = 0.4 ± 0.05. 

4.1 The cavity method 

The cavity method is a technique that allows one to calculate the average properties of 
statistical models denned on tree- like graphs in the thermodynamic limit j60j |6T]. The 
method, which is equivalent to the Bethe-Peierls (BP) approximation, is based on the 
assumption that due to the tree-like structure of the lattice, in absence of a given site (the 
cavity), the neighbors of that site are not correlated and their marginal joint probability 
factorizes. 

Removing one site from the graph creates C rooted trees, each one being a tree where 
all the sites of the bulk have the same connectivity C, apart from the root which has only 
C — 1 neighbors. The evaluation of physical observables is based on the determination 
of the properties of the site at the root of a rooted tree. Thanks to the above-mentioned 
factorization property one can write relatively simple recursion equations for the marginal 
probabilities of the rooted sites. Such equations have to be solved self-consistently, the 
fixed points of which yield the free energy of the system along with all the thermodynamic 
observables. 

The BP method can be interpreted either as the exact solution of the model on the 
tree, or as an approximate solution of the original model on the Euclidean lattice. In order 
to mimic an Euclidean lattice in D dimensions the tree should have connectivity C = 2D. 
Such an approach takes into account short-scale (0(1)) correlations and it is expected 
to give more accurate quantitative results than the standard fully-connected mean-field 
approach. It is well-known, for instance, that for the Ising model the BP method yields 
the mean-field critical exponents with a better estimate of T c than the one obtained on 
the fully-connected graph. 

4.2 The trees 

The definition of the vertices requires the selection of a particular orientation of the edges 
adjacent to a given site. We will define "horizontal" and "vertical" edges, each one with 
two possible orientations. With this procedure we associate a statistical weight to each 
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vertex configuration, even in such non-Euclidean geometry. In the recursion equations 
this partition will translate into four different species of rooted trees, depending on the 
"position" (left, right, down or up) of the missing edge. 

4.2.1 A tree of vertices 

In the models we are interested in, each site is a vertex and its coordination, which is fixed 
and equal to four, is the number of vertices connected to it. In order to distinguish one 
type of vertex from another, and to identify all possible phases, we define the analogue 
of the two orthogonal directions of the Euclidean square lattice: each vertex has four 
terminals that we call "up" (u), "down" (d), "left" (1) and "right" (r). So far, the vertices 
were labeled by their positions. Here, for the sake of clarity, we label them with a single 
latin index, say i, j, k, .... Vertices are connected through edges (i d j u ) and (i l k r ) that 
link respectively the down extremity of a vertex i with the up terminal of a neighboring 
vertex j, or the left end of i with the right end of its neighbor k. The symbols (i d j u ) 
and (i l k r ) denote undirected edges, so that (i d j u ) = (j u i d ) and (i l k r ) = (k r i l ). In this 
way, one creates a bipartition of the edges into horizontal (left-right (i l k r )) and vertical 
(up-down (i d j u )) edges. This notation gives a notion of which vertex is above (i) and 
which one is below (j) and similarly for the horizontal direction. See the sketch in Fig. [TUJ 
Each edge is occupied by an arrow shared by two vertices. An arrow defined on the 
(i l k r ) edge is the left arrow for the % vertex and the right arrow for the neighboring k 
vertex. A similar distinction holds for the vertical (i d j u ) edges. Each arrow, as any binary 
variable, can be identified with a spin degree of freedom, taking values in {—1, 1}. In this 
construction there are two kinds of spins, those living on horizontal edges, s/ji fe r\, and 
those sitting on vertical edg Without lose of generality, we choose a convention 

such that s/juid\ = +1 if the arrow points up and —1 otherwise. Similarly, s^ ri i^ = +1 
if the arrow points right and —1 otherwise, for the horizontal arrows (see Fig. [TUT) . This 
is the analog of the convention used for the spin sign assignment in the 2D model (of 
course, alternative choices of the signs of the spins can be used). 



S (i l k r ) — + 




Figure 10: (Color online.) Two spins living on an horizontal and a vertical edge, both taking 
value +1. 
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The local arrow configuration defines the state of the selected vertex. With the spin 
definition given above, we can assign a total spin Si to each vertex i, and define it as the 
sum of the spins attached to it, Si = \ J2j&di s {ij)> where di indicates the neighborhood 
of i. With this assignment, the spin associated to each type of vertices is Si = ±2 if the 
vertex is of type a, Si = if it is a b one, Si = if the vertex is of kind c, Si — for a d 
vertex, and Si = ±1 for the e ones. 
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Figure 11: (Color online.) Construction of un "up rooted tree" from the merging of a left, a 
right and an up rooted tree. 

Consider now a site (vertex) i at the root of a rooted tree, in absence of an edge 
with one of its neighbors, say site j. There are four distinct rooted trees depending on 
whether the missing edge with vertex j is the one on its left, right, up or down direction. 
By analogy with the 2D case, one could interpret these rooted trees as the result of the 
integration of a transfer matrix in four possible directions. This can be emphasized by 
taking into account the particular direction of the missing edge at the root that we will 
indicate as follows: i a — > j 13 , with a, j3 £ {u,d,l,r}. For instance, an "up rooted tree" is 
the one in which the root % has no connection to the up terminal of the vertex j, i.e. the 
link i d — > j u is absent. As shown in Fig. [TTJ such a rooted tree is obtained by merging 
a left, an up and a right rooted tree (with root vertices I, h and k respectively) with the 
addition of a new vertex % through the links V — > i l , h d — > i u , k l — > i r (pictorially the 
transfer matrix is moving down). Similarly, a "left rooted tree" is obtained by merging a 
down, a left and an up rooted tree, and so on. The Bethe lattice is finally recovered by 
joining an up, a left, a down and a right rooted tree with the insertion of the new vertex. 
Equivalently, given a tree, one creates rooted (cavity) trees by removing an edge. 

4.2.2 A tree of plaquettes 

As we will show in the following, the results obtained by using the tree structure described 
above compare extremely well to the ones in 2D in many respects. In order to further 
improve the approximation, in particular relatively to the nature of some transitions, we 
introduce a Bethe lattice of "plaquettes" (see Fig. [12]), where the basic unit cell is not a 
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single vertex, but a square of 2 x 2 vertices. This tree is constructed by connecting each 
plaquette to other four plaquettes (without forming loops of plaquettes) , in the same way 
as we constructed the tree of single vertices. In this procedure one has to be careful with 
the orientation of each plaquette and of the vertices on it, and with the order between 
the two edges outgoing from each side of the unit. 
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Figure 12: (Color online.) The main panel shows how to construct a Bethe lattice of individual 
units that can be chosen at will. In the right panel we show four different choices of such units. 
The image in (b) represents a single vertex that once inserted in (a) builds the simplest Bethe 
lattice of vertices, described in Sec. I4.2.T1 The tetrahedron in (c) is the individual unit used for 
the calculations of [62] and [53], as discussed in Sec. 14.2.31 Finally panel (d) shows a plaquette 
of four vertices that is the individual unit of the plaquette model, see Sec. 14.2.21 

In the following we will refer to the first simpler geometry as the "single vertex prob- 
lem" and to the second one as the "plaquette model" . 

4.2.3 Discussion 

A mean- field approximation for the pyrochlore spin-ice system in 3D, based on the same 
tree-like structure of single vertices described in Sec. 14.2.11 nas been already employed 
in [62] and [53]. In these papers no distinction between the orientation (up, down, left, 
right) of the edges was made. This approach was apt to deal with the SL and FM 
phases [53] only. Conversely, our approach keeps track of the four different directions and 
allows us to simultaneously study all possible phases, including the AF ones. Within our 
approach it is also easier to remove the degeneracy of the vertices by introducing external 
magnetic fields or other kinds of perturbations. 

The single-vertex BP approximation proposed in [62| 153] provides a very good quali- 
tative and quantitative description of the transition towards the frozen FM phase (KDP 
problem) in pyrochlore spin-ice in 3D. This is due to the absence of fluctuations in the 
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frozen FM phaseo However, such a single- vertex approximation is not precise enough to 
describe the unfrozen staggered AF order, due to thermal fluctuations. By constructing 
the tree of plaquettes we will manage to capture some of these fluctuations, and to obtain 
a more accurate description of the unfrozen phases. 

Let us finally mention that the ice-rule for the six-vertex model, or the "parity" rule 
for the eight-vertex might be viewed as a particular form of constraint that forces to 
flip all the variables on an entire loop in order to go from one allowed configuration to 
another. Such form of constraints has in many respects a strong algorithmic impact, 
in particular for algorithms that work with local updates. For this reason it has been 
investigated in the context of combinatorial optimization, for a wide class of problems 
mainly defined on tree-like graphs, with techniques similar to those adopted here for the 
tree of single vertices [63] . 

4.3 The six- and eight-vertex models on a tree of vertices 

In this Subsection we study the six- and the eight-vertex models on a tree of single 
vertices. We obtain the self-consistent recursive equations for the marginal probabilities 
along with their fixed points. We perform the stability analysis of the solutions, compute 
the free-energy of the different phases, and present the resulting phase diagram. 

4.3.1 Recursion equations 

We call ip^ 3 , with a,/3 G {u, d, I, r}, the probability that the root vertex i - in a 
rooted tree where (i a j^) is the missing edge - be of type ( G xl — { v ii ■ ■ ■ , v s}- Such 
probabilities must satisfy the normalization condition 

£ = 1 > v <* a /> ■ ( 16 ) 

In the recurrence procedure we will only be concerned with the state of the arrow on 
the missing edge. As a consequence, on the root vertex i of a rooted tree with a missing 
edge i a — > j 13 , we define ipf? = ip 1 ™^^ (+1) as being the probability that the arrow that 
lies on the missing edge (i a jP) takes the value s/ ia je\ = +1. Then, 

€ = + + + <^' u . 

# = €r jd + <^ jd + + <^ jd , 

4 = + + <r f + <~*' . ( i7 ) 

and 

= 1 - ^(+1) = 1 - 4 ■ (18) 

2 This is also due to the fact that D = 3 coincides with the upper critical dimension of the problem 
in the absence of defects (i.e., the six-vertex model): the mean-field BP approximation is almost exact 
in 31? apart from logarithmic corrections. 
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Clearly, other parameterizations are possible for these probabilities. For instance, one 
could use an effective field acting on the spin s^-m, i.e. = e h */(e h i + e~ h i), and the 
recursive equations would be equivalently written in terms of the fields h?. 

The operation of merging rooted trees allows us to obtain in a self-consistent fashion 
the set of probabilities associated to the new root in terms of those of the previous 
generation. In the bulk of the tree one does not expect such quantities to depend on the 
precise site, since all expected phases are homogeneous. As a result, the explicit reference 
to the particular vertex index % of the probabilities defined in eq. ( ITT)) can be dropped. 
By so doing, the following four coupled self- consistent equations for the probabilities ip a 
are obtained: 



tjj u = #«[ a> b, c, d, i> u , iff 1 , tp 
1 



g u (a,b, c,d,i/; u ,ip d ,ip l ,ip r )/z u 



+ 6(1 - - V r ) + c(l - ^")(1 - i) l W + dip l (l - i) u )(l - i) r ) 



Z" L 

ip l = &[a,b,c,d,ip u ,ip d ,ij; l ,ip r } = g l {a,b,c,d,ip u ,^ d ,ip l ,^ r )/z l 

{ ai) d ^ l i) u + 6(1 - - + c0 d (l - ^)(1 - i> u ) + d(l - i) d ){l - i) l )i) u 

iP d = #*»[ a , 6, c, d, ip u , ip d , ip r ] = 9 d (a, 6, c, d, ip u , ip d , V)/z° 

[ a^ r ip d ip l + 6(1 - i/> r )if> d (l - + c(l - ^ r )(! - i>V + # r (l - ^ d )(l - ft) 
^ = $r[ a , 6, c, d, ip u , ip d , V] = 9 r (a, 6, c, d, ip u , ip d , ip l , ip r )/z r 

[ ai> u i/; r i) d + 6(1 - i) u W{l - i) d ) + ci[) u (l - ^ r )(l - + d{l - 4> u )(l - 

(19) 

where a, b, c, and d are the fugacities of the vertices, as defined in Sec. [2] (see Figs. [TJ 
and E]), and z a are normalization constants that ensure the normalization of ip a : 



z a = g a (a, b, c, d, r, ^\ V) + 9 a (a, b, c,d,l- i[> u , 1-^,1-^,1- i[> r ) 



(20) 



The functions g a have been defined in Eq. ( TTTJ]) . The first term in this sum corresponds 
to the un-normalized contribution of a spin (+1) while the second term is for a spin (—1). 
For the sake of simplicity in eq. ( TT9|) we considered the argument in the functions ^> a to 
be the same for all the directions a, although in each equation the field defined along the 
opposite direction does not appear explicitly. 



4.3.2 Fixed points and free-energy 

In order to allow for a fixed point solution associated to c-AF and c?-AF staggered order 
we study eqs. (TT9l on a bipartite graph. We partition the graph into two distinct subsets 
of vertices A\ and A 2 , such that each vertex belonging to A\ is only connected to vertices 
belonging to A2 and vice versa. This amounts to doubling the fields degrees of freedom 
{i/>i, 1P2 h one f° r the sub- lattice A\ and the other one for the sub-lattice A 2 , and to solve 
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the following set of coupled equations: 



ij>2 = * a [a,b,c,d t tf,rl>t,rlt,rl>{] , 



(21) 



with a = u,d, r, I. The FM and PM phases are characterized by x/jf = while the AF 
phases by x/jf ^ ^ with if)" = 1 — . 

Considering the solution associated to A± only, the fixed points are 



V'PM = 


= (lp u = 




= hv- 


= = 


-\) 


VVFM 


= (lp u 


= 1,^ 


= i,r 


= i,^ 


= 1) 


^b-FM 


= (lp u 


= 1,^ 


= o,r 


= o,^ d 


= 1) 


V'c-AF 


= (lp u 


= 1,^ 




= o,^ d 


= 0) 


V'd-AF 


= (lp u 


= 1,^ 


= o,^ p 


= 1,^ 


= 0) 



(22) 



as can be checked by inserting these values into eqs. (TT9l) and (l2Tj) . The spin reversal 
symmetry implies that if}' = 1 — i/j is also a solution of the self-consistent equations. 
For the AF phase this corresponds to the excange of the two sublattices (i.e., the ex- 
change of ip\ with These solutions exist for any value of the vertex weights a,b,c 
and d. Conversely, the stability of the solutions, that will be discussed in more detail 
in Sec. 14.3.41 depends on the fugacities. The numerical iteration of the self-consistent 
equations, eqs. ffT9|) and fl2T|) . confirms that the fixed points given in eq. fl22|) are the only 
possible attractive solutions in the different regions of the phase diagram. 

In order to calculate the free-energy, it is useful to consider the contributions to the 
partition function coming from a vertex, an horizontal edge and a vertical edge. These 
quantities are defined as follows: 



a 

+ b 
+ c 



^ u ^ r ijj d + (1 - i) u )(i - i) l )(i - i) r ){\ - i) d 
(1 - ip l )ip u (i - 4) r )ip d + i/> l (i - ip u )^ r (i - ii) d ) 
(1 - i) u )(i - 4) l )^ r %i) d + ^V(i - ^ r )(! - i> d ) 

d - ip u )(l - 4j r )tp d + (1 - x/j l )^ u ^ r (l - tp d ) 



(23) 



and 



V } [^,^]=^ + (1-^)(1 



(24) 



The first term Z v represents the shift in the partition function brought in by the in- 
troduction of a new vertex which is connected with four rooted trees. The other terms 
Zn r \ and Z/ U( a represent the shift in the partition function induced by the connection of 
two rooted trees (respectively one left and one right or one up and one down) through 
a link. In terms of these quantities one can compute the intensive free-energy (here and 
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in the following we normalize / by the number of vertices) which characterizes the bulk 
properties of the system in the thermodynamic limit: 



(3f[a,b,c,d, iff 1} if> 2 ] = --(lnZ^x] +\nZ v [i/> 2 



- In Z {lr) [^U 2 r ] - In Z {lr) [if, l 2 , ^] - In Z {ud) [tf, ^} - In Z {ud) fo$, 

The free-energy (125]) evaluated at the fixed points (122]) reads as follows: 

a + b + c + d 



(25) 



(26) 



Pfpu = Pf[a, b, c, d, V> PM ] = - In 

PfarFM = Pf[d, 6, C, d, VVfm] = ~ ln a 

/3/fe-FM = /3/[a, b, c, d, V'^fm] = - ln& , 

/3/c-af = Pf[a, b, c, d, V c -af] = - In c , 

0fdrAF = Pf[a, b, c, d, V^-af] = - In d . 

By comparing these free-energy densities one readily determines the phase diagram 
along with the critical hyper-planes separating the different phases. For instance, from 
the condition f PM = / n _ FM one finds a c = b + c + d. Surprisingly enough, it turns out 
that the critical planes have exactly the same parameter dependence as in the six- and 
in the eight-vertex model on the square lattice. The BP approximation yields, therefore, 
the exact topology of the phase diagram. 

4.3.3 Order parameters 

According to the definitions in eq. (JT]), we characterize the phases by direct and stag- 
gered magnetizations. In particular, we define the following order parameters, each one 
associated with a particular ordered phase and all of them vanishing in the PM stateS 



Z 



V 



m b-FM = i(l - ^)^(1 - W ~ ^(1 - W(l - V)\ 



J v 
C 

Z v 
d 



m c _ AF = — [^V(i - ^ r )(i - i) d ) - (i - ^)(i - 



(27) 



m d .AF = ^r[(l~ ^)^> r (l - V) ~ ^(1 - ^)(1 - VW] , 

Zj v 

where Z v is the contribution of one vertex to the partition function defined in eq. 
The order parameters 027p are easily evaluated at the fixed points, eq. f[2"2"]) . The PM so- 
lution ifjpM yields vanishing magnetizations for all the order parameters. Conversely, any 



3 This is a generalization of the definitions given in eq. (TTJ), needed in order to make the difference 
between FM orders dominated by a or b vertices, as well as AF orders dominated by c or d. It corresponds 
to consider m a _FM = \{ m °+ + m +)i m b-FM = \{ m + — to +)j m c -FM = \ijn x _ — m v _) and md-FM = 
^{rrfL + m v _). 
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ordered solution i/'fp=* gives a saturated value of the associated magnetization mf p= * = 1 
(and zero for the other components mf p ^* = 0), corresponding to a completely frozen 
order phase. 



4.3.4 Stability analysis 

The stability of the fixed points is determined by the eigenvalues of the Jacobian matrix 



M 



d^ 



, , (28) 



that describes the derivative of the vector function = (^ u , defined in eq. f[T9"|) 

with respect to the fields {ijj a }, evaluated at the fixed point, t/> fp . The eigenvectors are 
parametrized as (Si/j u , 5ip r , 8ip l , 8ip d ) . The solution if) { becomes unstable as soon as one 
of the eigenvalues becomes one, |^ maa; | = 1. 

Indeed, one can easily prove that the condition |£? max | = 1 corresponds to the di- 
vergence of the magnetic susceptibility and allows us to identify the different phase 
transitions. More precisely, the susceptibility reads 

X = N dh = N I* {s ^ s ^ )c x 2^ 2^ 2^W) C (29) 

(ij),(kl) a,l3=(ud),(lr) r=l V{r) 

where in the last equality we used the homogeneity of the solution. The symbol V(r) 
indicates that the sum runs over all the paths that connect a given spin on an edge of 
type a, supposed to be the centre (site denoted by 0) of the tree, to all the spins that 
live on edges of type /3 and located at a distance r from (in terms of the number of 
edges that make the path). As the tree has no loops, such paths are uniquely defined. 
The above formula can be simplified by first using the fluctuation-dissipation relation 

where Hq is a field conjugated to Sq, and next using the chain rule 

^ = &^ r r ^d^ L d(4l (si) 
dh% d^o f_| d^ 71 " 1 dtpyr-i ' 1 ' 

where the particular values taken by {7j} G {u, d, I, r} depend on the path followed. Each 
derivative is finally evaluated at the fixed point. Then, defining the vectors \v a ) such that 
v a = ffs - ' an d \ w a) such that w@ = ^r, one obtains 

oo 

a,P=(ud),(lr) r=l V(r) 

^ Z^Z^ d^ -1-1 d^-i d^- 1 ^ 

a,/3=(ud),{lr) r=l V{r) ru i=2 r r 
oo oo 

oc ( v a\M r - 2 \w^) ~ ^TrilT. 

a,/3 = (ud),(lr) r=\ r=\ 
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As long as the absolute value of the eigenvalues remains smaller than one, i. e. \E„ 



< 1, 



the series converges yielding a finite susceptibility. Otherwise, if \E„ 
diverges yielding an infinite susceptibility. 



> 1 the series 



Let us focus on the stability of the PM phase V'pm' ? - e - ^ 
matrix M evaluated in the PM solution is 



The 



M, 



d 



d 



b + c — d 



a+b+c+d a+b+c+d a+b+c+d 
a—b+c—d a+b—c—d 



a+b+c+d 

—a + b + c — d 
a+b+c+d 



a+b+c+d 



a+b—c—d 







-a + b + c — d 
a + b + c + d 



a + b + c 
a — b + c 



d 
d 



—a + b + c — d 
a+b+c+d 

a—b+c—d 
a+b+c+d 

a+b—c—d 
a+b+c+d . 



(33) 



Its eigenvalues are 

E PM 



3a — b — c — d 



E. 



PM 



a 

a - 



-b + c- 
b- 3c 



d 

-d 



E. 



PM 



-a + 3b 



c 



b + c + d 



E 



PM 



a+b+c+d 
a + b + c — 3d 



(34) 



a 



c 



d 



In the order (Sip u , 5ip r , Sijj 1 , 5ijj d ) the corresponding eigenvectors can be written: Vi = 
(1, 1, 1, 1), v 2 = (1, -1, -1, 1), v 3 = (-1, 1, -1, 1) and v 4 = (-1, -1, 1, 1). In general, the 
eigenvalue Ef M regulates the instability towards the a-FM, E% M towards the fe-FM, Ef M 
towards the c-AF, and Ej M towards the d-AF. 

One reckons that the stability of PM solution can be stated in terms of the condition 



Eg M 



)(1 



E™) 



[1 -Ef M ){l 



E™] 



E™] 



[1 + E™)(1 
which in terms of a, b, c, d reads 

a 2 + b 2 - c 2 - d 2 



2(ab + cd) 



E™){l-E™) 



|A 8 | < 1 . 



< 1, 



(35) 



(36) 



Therefore, the stability analysis of the PM phase also yields the same condition for the 
phase boundaries as in the exact solution of the 2D model on the square lattice [see 
eq. © and Fig. [2]. 

A similar analysis can be carried out to evaluate the stability of the FM and AF phases. 
One has to evaluate the matrix at the ordered solutions' fixed points. Considering 
the fixed point i/Vfm as an example, one obtains the four eigenvalues: 



El 



a-FM 



E. 



a-FM 



b + c + d 
a 

b — c + d 
a 



E. 



a-FM 



C 



d 



E 



a-FM 



a 

c — d 
a 



(37) 
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and the corresponding eigenvectors are the same (in the same order) that we have already 
discussed for the PM solution. From eq. (1371) . one sees that E®~ FM = 1 as soon as 
a = b + c + d, i.e. when A§ = 1, while for larger values of a the a-FM fixed point is 
stable. Moreover, the solution may develop multiple instabilities when some vertices are 
missing. Analogous results hold for the other ordered phases. 

Finally, let us remark that the model also has discontinuous transition points sepa- 
rating ordered phases, for instance the a-FM and the c-AF at a = c and b = d = 0, 
the 6-FM and the c-AF at b = c and a = d = 0, the c-AF and the d-AF at c = d and 
a = b = 0, etc. 

4.3.5 The six-vertex model: phase diagram and discussion 

We found four different fixed points in the six- vertex model (d = 0). These characterize 
the (i) a-FM phase, (ii) 6-FM phase, (iii) c-AF phase and (iv) PM phase (we will distin- 
guish between the PM phase found on the tree and the SL phase in D = 2 in ways that 
we will describe below). 

The transition lines were found in two equivalent ways. On the one hand, we compared 
the free energies of the PM and ordered phases as given in eqs. f[2"5"j) . On the other hand, 
we analyzed when one of the eigenvalues of the stability matrix becomes equal to 1, 
signaling the instability of the considered phase. We found that, as in the exact solution 
of the model in 2D, the phase transitions are controlled by the anisotropy parameter 
A 6 defined in eq. (J4]), with the transition lines determined by |A 6 | = 1. Interestingly 
enough, the eigenvalue E\ M equals 1 when d = 0, V a, b, c, meaning that the whole 
PM phase is critical in the sense that the magnetic susceptibility diverges for d = 0. 
This is reminiscent of the critical properties of the SL (or Coulomb) phase in D = 2. 
Analogously, if a = then Ef M = -1 V c, 6, d; if b = then E% M = -1 V a, c, d and 
if c = then Ef M = 1 V a, 6, d. These observations suggest that in fact the PM phase 
becomes critical whenever one of the four kinds of vertex is absent. 

The transitions between the PM and the ordered phases are all discontinuous (as 
it can be seen from the singular behavior of the free-energy): the (possibly staggered) 
magnetization jumps from zero to one. The fact that the magnetization saturates at its 
maximum value in the whole ordered phases indicates that in the FM and AF phases the 
order is prefect and frozen. 

Even though the transitions are discontinuous, they are characterized by the absence 
of metastability and hysteresis, and they are associated to a diverging susceptibility. In 
fact, approaching the transition from the two sides, both the PM and the FM solutions 
become unstable. This kind of transition corresponds to a multi-critical point of infinite 
order and it is called KDP transition in the literature [51] . Note that if one focuses on the 
transition line, and plugs the critical value a c = b + c into eqs. ffTU]) and (|2"T1) . assuming 
the homogeneity of the solution, ip a = ip Va, the equations take the trivial form ip = if), 
meaning that all values of ip are solutions of the self-consistent equations and all values of 
magnetization are equally probable. Actually, the free-energy for the same critical value 
of a c does not depend on ip and it is thus the same for all values of the magnetization. 
The PM-AF transition shares the same properties of the KDP transitions between the 
FM and the PM phases. 

At the spin ice point, a = b = c, the entropy is S sv /N = —f3f sv = In 3/2 ~ 0.405, i.e. 
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the Pauling result for water ice j9] (see Fig. [T6|) . 

4.3.6 The eight-vertex model: phase diagram and discussion 

The addition of d vertices does not change the properties of the a-FM, b-FM and c-AF 
phases. As vertices of kind d are allowed, another AF phase, denoted d-AF, emerges, 
corresponding to a phase with staggered order of four-in and four-out vertices. The 
transition from the PM to the d-AF phase is also discontinuous and the ci-AF phase is 
completely frozen as well. Moreover, as for the six- vertex model both the FM and the PM 
solutions become unstable on the transition planes. When all the fugacities are finite the 
four eigenvalues of the stability matrix within the PM phase become all smaller than one 
(except, of course, on the critical transition planes). This means that the PM phase of the 
eight-vertex model is no longer critical. The transition lines are altogether characterized 
by the condition |Ag| = 1, as for the exact solution of the 2D model, with A§ given in 
eq. 09]). 

4.3.7 Summary 

In short, the location of the transition lines of the 2D six- and eight-vertex models are 
reproduced exactly by the calculations on the tree of single vertices. However, the critical 
properties on the tree are different from the ones of the 2D model. The absence of loops 
'freezes' the ordered phases and makes all transitions discontinuous, with the (possibly 
staggered) magnetization equal to one in the whole ordered phases. While this is true for 
the PM-FM transition of the six- vertex model in 2D, this is no longer true neither for 
the PM-FM transition of the eigth-vertex model, nor for the PM-AF transitions of the 
six- and eight-vertex models (in particular, in the former case the PM-AF transition is 
of KT type). 

4.4 The six- and eight-vertex models on a tree of plaquettes 

In this Subsection we study the six- and eight- vertex models on a tree of plaquettes. The 
calculations proceed along the same lines as for the single vertex tree. They become, 
though, more involved, since the number of configurations allowed on a plaquette is quite 
large. 

4.4.1 Recursion equations 

Each rooted tree of plaquettes has two missing edges, which means that one has to write 
appropriate self-consistent equations for the joint probability of the two arrows lying on 
those edges. The analogue of ip a , which in the previous section described the marginal 
probability of the arrow to point up (or right, depending on whether we are considering a 
vertical or an horizontal edge), now becomes a probability vector with four components. 

The marginal probability to find a pair of arrows with values ++, — h, H — , will 

be denoted by V° = {>++, V>--}0 

The superscript a = u,d,l,r denotes, just 

4 In the case of the tree of simple vertices the bold symbol referred to the vector associated to the 
four directions a, while here it is signaling that already for one direction one has to deal with multiple 
probabilities. 
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as before, whether the pair of arrows are on the missing edges of an up, down, left, 
right rooted tree. This is illustrated in Fig. [131 We keep the convention on the choice 
of the signs of the spins; namely, the spins take positive values if the arrows point up 
or right. Moreover, we assume that the first bottom index indicates the state of the 
arrow that is on the left, for the vertical edges, and on the top, for the horizontal ones. 
Consequently, the second bottom index refers to the value taken by the spin sitting on 
the right (respectively the bottom) edge for vertical (respectively horizontal) edges. 
The weights of the vertex can be written as follows 

w S i,s2,s 3 ,s 4 ( a ,b,c,d) = - a'{l + s 1 s 2 s 3 s< l ) + b'(s 1 s 3 + S2S 4 ) + c'{s 1 S4 + S2S 3 ) + d'(s 1 S2 + s 3 s 4 ) 

(3* 

where s\, . . . , s 4 are taken as in Fig. HU and 

a' = -(a + b + c + d) , b' = -(a + b — c — d) , 



2 y ' 1 2 

d = - (a — b + c — d) , d' = - (a — b — c + d) 



(39) 



u,d 



l.r 



Figure 13: Definitions of ^--ja^^.r.d used in the recursion equations for the 

plaquette model. First line: ij) u,d where the first index ± denotes the value of the spin on the 
left and the second index denotes the spin on the right. Second line: tp ' where the first index 
± denotes the value of the spin on the top and the second index denotes the spin on the bottom. 

Similarly, we introduce a new parameterization of the probability vector, 



a 



[l + s iS 



[Si + s, 



+ (1 - SiSj 



+ (Si - s, 



2 



1 

4 L 



1 + SiSj s a + (si + Sj) p a + (f 
where we introduced a new set of variables defined as 



4> a = ( P a , s a , q a ) = - + c_ - c + - ^ 



(40) 



(41) 



We exploited the fact that, due to the normalization conditions, only three variables 
are truly independent for each edge direction. In the PM phase there is no symmetry 
breaking and p a = q a = 0. The ordered phases and the corresponding phase transitions 
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Figure 14: Left panel: representation of the variables {sj} used in the definition of the vertex 
weight in eq. (|38|) . Right panel: numbering assigned to the spin/arrow variables in eqs. (|43p 
and ([4"5]) . We denote by Sp = {s\, S2, S3, S4, S5, sq, S7, ss, £i, £2, £4} the set of spin variables on a 
plaquette. 



are characterized by the spontaneous symmetry breaking associated to FM, p a 7^ 0, or 
AF, q a 7^ 0, order respectively. 

For the sake of completeness we also report the inverse mapping: 



1 

4' 
1 

4< 



= - 7 {l + s a + 2p c 



= -y-(l-s a -2q a ) 



(42) 



#L = -(l + s a -2p° 



The self-consistent equations for the probability vector describing up, left, right and down 
rooted trees now read 



OC W ^^2,t2A^2, S 3, S 4,t3 W t4,t3, S5 , S6 ^ 8 A ! t4, S 7'0UlCr S6 '0r i 



5S4 ' 



W Sl,S 2 ,t2,tl' lt 't2,S3,S4,t3' lt '*4,t3,S5,'S6 

Sp\{s5,SA} 



S 7 S 6 5 



(43) 



OC ^ ^ 1 , S 2,t2A W t2, S 3, S 4,t3 W t4,t3, S5 , S 6^ 8 A ! t4, S 7'0UlV'f 2S 3^l 

Sp\{s 6 ,s 7 } 

OC ^ ^l,S2,<2,il W i2,S3,S4,i3 U; i4,i3,S5,S6 U; S8,tl,i4,S7 ? / ; r 7 S6^5S4^S2S3 ' 
Sf>\{ s 8,si} 



S5S4 ' 
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and the normalization constant is given by z a = ^2 S . s . i/)f. Sj - 

It is more convenient to work with the variables <fi a = ^ , 0f ) = {p a , s a , q a ) for 
which one can readily derive the set of self-consistent equations from eqs. fj4"3l) . 



p a = <£?[a, 6, c, d, <f>\ 0', r , </> d ] = n + - , 

s a = $«[ a , 6, c, d, U , 0', r , d ] = + - - if 



+- ' 



(44) 



with a = u,l,r,d. The argument of the functions of the right hand side can be 
expressed in terms of the variables p a , s a , and q a using the transformations 

4.4.2 Free-energy 

The free-energy per vertex can be generically written as: 

(3f[a, b, c, d, r, <4> r , ¥] = \ (in KMn + ln E ~ ln Z vi ) 



Sl,S 2 



with the plaquette partition function given by 



Zpl — E Ws l' s 2,t2,il^i2,S3,S4,i3^4,t3,S5,S6^ 8 ,h,i4,S7V ; s 8Sl ^r 7 S 6 ' ? / ; l5 S4 V ; , 



«2«3 ' 



(46) 



where one can rewrite ?/> Q in terms in (f> a via eq. 



4.4.3 Order parameters 



The definitions of the order parameters fl27j) given in Sec. 14.3.31 have to be generalized 
to take into account the Sp variables. The order parameter associated to the a-FM 
transition reads: 

«vfmK b i °i d i i> u , *l>\ V, if> d ] = 

^pl ^ - 



S P 



(47) 

with O^fmISp] = f (X)i=i,...,8 f + Si=i,...,4*iJ an d ^ defined in eq. (T46]). The other 
order parameters are obtained from eq. (I4"?j) by substituting the quantity O a -FM with other 
staggered magnetizations suited to capture the different ordered phases. More explicitly, 
they are 



Ob-FM[Sp] 

C c -fm [Sp] 

Od-FM[Sp] 



L 2 ( _Sl + 82 + s 3 - s 4 - s 5 + s 6 + s 7 - s 8 ) + *i - *2 + h - U 



1 

L2 
1 

2 



— Si + S2 — S3 — S4 + S5 — Sg + S7 + Ss) — ti + ^2 + ^3 — ti 

Si + s 2 - s 3 + s A - s 5 - s 6 + s 7 - s s ) -t 1 -t 2 + t 3 + t 4 



(4* 
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4.4.4 Stability analysis 



Similarly to what done for the model on the tree of vertices, we investigate the stability 
of the solutions by studying the stability matrix 



M 



a,/3 



d$? 



i -j 



a, (3 = u,l,d,r , 



i,j = 1,2,3 



(49) 



In general, this is a 12x12 matrix but, depending on which solution one is studying, 
it might break into disjoint blocks. Moreover, the entries of such matrix have certain 
symmetry properties that simplify the calculation. From the stability analysis it turns 
out that the values of the fugacities where the PM solution becomes unstable also coincide 
with the transition planes between the different phases, corresponding to the points where 
the free-energies of the different phases cross. We adopt the following parametrization of 
the 12 components of the eigenvectors 



6<f = (5p u , 6s u , 5q u ; 5p d , 5s d , 6q d ; 5p l , 5s l , 5q l ; 5p r , 6s r , 5q r 



(50) 



4.4.5 Fixed points 

Let us discuss each fixed-point solution - phase - separately. 

The paramagnetic phase. The PM fixed point is of the form PM = 4>p M = 4>vm = 
0pm = 0pm = (Ppm = 0, spm, 9pm = 0) with spm determined by eq. (jSJ). For convenience 
we introduce the quantities 

a — b 



c — d 



x 



y 



a+b+c+d a+b+c+d 

These relations can be easily inverted as 

a 1 , 

-(l + z + 2x) , 



a + b — c — d 
a+b+c+d 



(51) 



a+b+c+d 

c 

a+b+c+d 



z + 2y) 



b + c + d 
_d 

b + c + d 



4 



{1 + z- 2x) 



id 



(52) 



2y) 



Thus the fugacities of the vertices are linear functions of x, y and z up to a normalization 
factor (a + b + c + d)^ 1 . The self-consistent equation for spm now takes the form 



{a + b + c + d)\x 2 -y 2 )(l + z 2 ) 
with 



T + 1 



T + 1 



SPM 



b PM 



T(a, b, c, d) 



T-l riV1 T-l 
Ax 2 + z 2 -l (a + b)(c + d)-(a-b) 2 



s 4 

5 PM 







Ay 2 + z 2 -l (a + b)(c + d) - (c - d) 2 ' 
Apart from the trivial solutions spm = 1, — 1, the PM solution is given by 



■SPM 



T 



(53) 



(54) 



(55) 
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Equation (J55J) implies that the PM solution depends upon a single parameter T, and that 
it remains unchanged if the vertex weights a, b, c, d are modified without changing the 
value of T. The limit of infinite temperature of the eight- vertex model, i.e. a = b = c = d, 
corresponds to spm = which implies = if)"_ = ip+_ = ip® + = 1/4, Wa G {u, I, r, d}. 
We anticipate that this result will also be found in the sixteen-vertex model. On the 
contrary, the limit of infinite temperature of the six-vertex model a = b = c and d = 
(more generally, whenever one of the four kinds of vertices is absent) corresponds to a 
non-trivial solution spu 7^ which implies = ^ = 

By inserting the PM solution in the free-energy density (|4"5]1 one obtains 



/3 f PM (a,b,c,d) = -In 



a + b + c + d 
2 



1 r (x 2 -y 2 ) 2 (-3 + 2x 2 + 2y 2 - z 2 ) 
In 



4 



2x 2 + 2y 2 + z 2 



(56) 



where x, y and z are given in terms of the vertex weights in eq. (I5ip . This function is 
clearly invariant under the exchange of a with b, c with d, and the simultaneous exchange 
of the weights of the FM vertices with the AF ones. 

Ferromagnetic phases. In order to study the ordered phases we introduce the following 
functions of four variables: 



wl-w 2 - (w 3 + W A 

f % = — 2 7 \2 ' 57 

w{-w{ - (w 3 - w 4 ) 2 



(59) 



X(I W \) = i\n\ W.-W.-W.-W, 

11 m]) 4 I2w 2 w 2 (w 2 + w 2 ) + K + w 2 w 2 ){w 2 - w 2 - w 2 - wf 
with {w m } = W\, W2, W3, W4, and 

(wl-w%-(w 3 - w 4 ) 2 ) (wf -vol - {w* + W4) 2 ) 

H\{w m \) = ======== 

(wj — w\ — w\ — w 2 ) \J (wf — w\ — w\ — w\) 2 — Aw 2 w 2 

(2w1w\ + wf (w 2 — wj—w 2 — w\))w 2 
(wf(wf(wf — w\ — w\ — w 2 ) + 2w 2 wj) — w 2 w 2 (—w 2 — w\ + w\ + w 2 )) 

Note that z/, E and n are symmetric under the exchange of w 3 and W4. 

The a-FM phase is homogeneous along all the directions and it is given by 

0a-FM = 0a-FM = (pa-FM = \J v{d,b\ C, d), S a . FM = 1, g a _ FM = Oj , (60) 

Va = u,l,r,d. This implies that = ip |_ = 0. Spin reversal symmetry is sponta- 
neously broken since ^ ip The associated free-energy and magnetization read 

P /a-FM (a, b, c, d) = E(a, b; c, d) , 

(61) 

m a .FM(fl, b, c, d) = /i(a, b; c, d) . 

The extension of these results to the other FM phase is straightforward. The b- 
FM is characterized by the same functions, with the exchange of b and a. One finds 
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^b-FM = \Pb-FM = —\/v(b, ®] c, d) , s 6 _fm = 1, Qb-FM — Oj for the horizontal edges while 
the positive sign remains on the vertical ones. 

Antiferromagnetic phases. AF order can also be characterized by the functions u, £ 
and /i defined in eqs. fl57|) . fl58l) and fl59l) . In particular, for the c-AF phase, 

0c-AF = 0c-AF = (Pc-AF = 0, S C -AF = ~1, <?c-AF = \/l>{c,d\ (2, 6) J . (62) 

This implies ^ ++ = -0 =0 and a staggered order with + _ 7^ "0-+- Spin reversal 

symmetry and translational invariance are simultaneously broken. The free-energy and 
staggered magnetization read 

P /c-af (a, b, c, d) = S(c, d; a, 6) , 

(63) 

m c _AF(a, c, cZ) = /i(c, d; a, 6) . 

The AF solution dominated by d vertices is of the same form as the one above, with the 
exchange of c and d, and q^^p = —Qc-af meaning that the c-AF and d-AF solutions only 
differ by a sign along the horizontal edges. This is due to the fact that a d vertex can be 
obtained from a c vertex by reversing its horizontal arrows. 

4.4.6 The six-vertex model: phase diagram and discussion 

The PM solution, eq. (|55|) . evaluated at d = describes the PM phase of the six- vertex 
model. The fixed point is characterized by a value of spm in the interval — 1 < spm < 0. 
As for the model on the tree of vertices, as soon as one of the fugacities is set to zero, 
the stability matrix (I49p has an eigenvalue whose absolute value is equal to one in the 
whole phase. This calculation allows to show explicitly how the introduction of a hard 
constraint can drastically change the collective behavior of the system: in our case, the 
addition of a hard constraint turns the conventional PM phase into one with a diverging 
susceptibility and a soft mode. The normalized eigenvector associated to this mode is of 
the form 

Sippu = (1/2, 0, 0; -1/2, 0, 0; -1/2, 0, 0; 1/2, 0, 0) . (64) 

The PM-FM transitions are of the same type as the ones found in the single vertex 
problem. Eqs. flB"Ul) and (IBTl) evaluated at d = yield a discontinuous transition towards 
a frozen phase with p a -FM = 1, i-e- ^++ = 1 5 at a c — b + c (or b c — a + c) where spm = 
V6, c, accordingly to eq. (155|) . The free-energy density in the frozen phase, f a -FU = — In a, 
and the magnetization, m a _FM = 1, are identical to the exact results on the square lattice. 
At the transition, /3/pm(q c , b, c) = /3/ a -FM(a c , b, c) = — lna c . 

On the transition lines both the PM and FM solutions become unstable. One can 
check that fixing a = b + c (or equivalently b = a + c) and looking for a solution of the 
type 0a_ fm = (P'a—FMi s Lfm- ?«-fm = °) the two equations for p c a _ FU and s c a _ FM become 
linearly dependent and simultaneously satisfied by the condition (Pa_ FM ) 2 — s a _ FM = 0, 
that describes an entire line of fixed points joining the FM solution (p a — fm = 1, Vfm = 
1) to the PM one Q»pm = 0, spm = 0). In Fig. [15] we show the free energy in the plane of 
the solutions (s,p), with s a = s, p a = p and q a = Va, and the free energy diminishes 
from light to dark colors. The three panels correspond to different values of the fugacities 
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and represent from left to right a paramagnetic minimum \Aq\ < 1, the a-FM critical 
point A6 = 1 associated to a degenerate line of minima, and the a-FM phase with Aq > 1. 




Figure 15: (Color online.) Projection of the free energy in the plane (s,p, q = 0), where p a = p, 
s a = s and q a = Va, for some value of the fugacities for the six-vertex model. The free energy 
decreases from light to dark colors and the minima are in the black regions. The three panels 
show from left to right the PM phase a < b + c, the critical point a = b + c and the FM phase 
a > b + c. The figure shows that at the transition point the free energy has an entire line of 
degenerate minima. 



The PM-AF transition for the tree of plaquettes is still placed at c c = a + b (as for the 
tree of single vartex and for the exact result) but it is qualitative different from the one 
found with the single vertex. The small loops of four spins allow for fluctuations in the 
AF phase and the transition becomes a continuous one with a singularity of the second 
derivative of the free-energy. One can indeed check that /pm(O) b, c c , 0) and /c-af(o, b, c c , 0) 
have the same first derivatives with respect to the fugacities at the critical point. At the 
AF-PM transition the PM solution reaches the critical value spm = — 1 and for larger 
values of c the PM solution becomes imaginary. The magnetization, m c _AF, is given by 



V2(a + bf[(a + b) 2 + ab] lc 

m c _AF(c; a, b) = —= * h O 

V^b[{a + by + ab({a + b) 2 + a&)] V 



C - C c ,3/2 



(— ) 



(65) 



close to the transition line, which gives the mean-field classical exponent (3 = 1/2 (see 
Fig. ED . 

In Fig. [16]we compare the free-energy of the model on the tree of single vertices (f sv ), 
on the tree of plaquettes (f p i), and the exact free energy of the model on the 2D square 
lattice (/21)) pT]- The left panel of Fig. [M] shows the free-energy in the PM and AF 
phases as a function of a/c, moving along the line a = b. The figure clearly shows that 
the discontinuity of f sv at the transition is smoothed out by the inclusion of small loop 
fluctuations. In the right panel we show the free-energy in the PM and FM phases, as a 
function of a/c moving along the line ab = c 2 . 

In the spin-ice point a = b = c the entropy per vertex of the six-vertex model on a 
tree of plaquettes is 

S*/N = -(3f pl = ~ In A ~ 0.418 . (66) 

This result is closer to the exact value obtained by Lieb for the 2D model S 2D = 
(3/2) ln(4/3) ~ 0.4315 (see Fig. QSJ. 



38 



As shown in Fig.[16]the free-energy of the frozen FM phases is the same for the square 
lattice model, the tree of single vertices and the tree of plaquettes. Instead, in the PM 
and AF phases the mean- field approach does not give the exact result. The plaquette 
geometry yields a better quantitative estimation of the thermodynamics properties of the 
model compared to the single vertex geometry (f 2D < f P i < /s<u)-@ 

The properties discussed so far are a general consequence of the form of the functions 
defined in eqs. ( 157|) . (|58|) and ( |59|) . independently of the precise specification of their 
arguments. As a consequence, similar conclusions are drawn when any of the four types 
of vertices is missing. Note that for yj\ > w 2 , w 3 and w<± = one recovers: 

u(w 1 ,w 2 ;w 3 ,0) = h>(w 1 ,w 2 ;0,w 4 ) = 1 , 
E(u>i, w 2 ; w 3 , 0) = £(iwi, w 2 ; 0, iu 4 ) = - In w\ , 
fj,(w 1 ,w 2 ;w 3 ,0) = 1 , 

while for w 2 = and W\ > w 3 , W4 ^ they all take a non trivial value. This means 
that whenever one of the fugacities of the AF vertices (c or d) vanishes, the FM phases 
are frozen (since /i = 1 and £ = cte) while the AF transition is continuous. The same 
holds for the AF order, which becomes frozen, in absence of one of the FM vertices [see 
eq. dSI])]. 




Figure 16: (Color online.) Free-energy density of the six-vertex model along the two paths 
sketched in the insets (red dotted lines), a/c = b/c (left) and ab = c 2 (right). Left panel: 
free-energy in the PM phase (a/c > 1/2) and AF phase (a/c < 1/2) on the Bethe lattice of 
single vertices (red curve) and plaquettes (green curve), and of the exact results on the 2D 
model (pink and blue lines). Right panel: free-energy of the PM phase for 0.6 < a/c < 1.6 on 
the Bethe lattice of single vertices (red curve) , of plaquettes (blue curve) and in the 2D model 
(pink line). In the FM phases, both single vertex and plaquette trees lead to the same (exact) 
free-energy (green and light blue lines). 



5 The mean- field approximation can be interpreted as a variational principle on the free energy of the 
system. As a result, a more accurate mean-field approximation yields a smaller free energy, closer to the 
true one of the 2D system. 
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a/c a/c 

Figure 17: (Color online.) The staggered magnetization, m c _AF, on a tree of single vertices 
(green curves) and of plaquettes (red curves). Left panel: the six- vertex model with a = b. 
The green step shows the transition towards a completely AF ordered phase (single vertices) 
while the red curve shows a continuous transition (plaquettes), see eq. ()65|) . The inset shows 
the mean-field exponent /3 = 1/2, m^_ FM ~ [(a — a c )/c]. Right panel: the eight-vertex model 
with a = b and d/c = 0.2. PM-to-frozen AF transition on the single-vertex tree (green) and 
discontinuous behavior on the plaquette tree (red), see eq. (follj) . showing a jump towards the 
AF phase (with a non- frozen staggered order due to thermal fluctuations). The dashed red 
line shows the continuation of the FM solution (|6ip beyond the transition point where the PM 
phase is favored. 



4.4.7 The eight-vertex model: phase diagram and discussion 

For the eight-vertex model the plaquette geometry gives quite different results with re- 
spect to the tree of single vertices. Indeed, on the tree of single vertices the addition of 
vertices of type d does not change the nature of the transitions (which are all discon- 
tinuous and of critical-to-frozen type as explained in Sec. 14. 3p . whereas on the tree of 
plaquettes it does. However, even though using the plaquette geometry the nature of the 
phase transitions are changed with respect to the tree of single vertices, their location is 
not modified (the location of the critical planes will still coincide with the exact solution). 
On the other hand, as far as the PM phase is concerned, both within the single vertex 
mean-field approximation and the plaquette one, the criticality and the "soft modes" 
disappear as soon as a, b,c,d ^ 0. 

When a,b,c,d ^ the free-energy at the transition planes shows a singularity in its 
first derivatives corresponding to a first-order phase transition. Indeed, one can check 
that 

fpM(w 2 + w 3 + w 4 ,w 2 ;w 3 ,w 4 ) = T,(w 2 + w 3 + w 4 ,w 2 ;w 3 ,w i ) , (67) 

where w 3 and are the statistical weights of FM (resp., AF) vertices if the transition 
under consideration is a PM-FM (resp., a PM-AF) one. Consequently, the magnetization 
at the transition shows a finite jump towards a non- frozen ordered phase. The a/c 
dependence of the staggered magnetization m c _AF is displayed in the right panel of Fig. [T7] 
for the tree of single vertices (green curve) and plaquettes (red curve). 

Let us now analyze the a-FM-PM transition and focus on the transition plane, a c = 
b + c + d. If we plug a solution of the type 0q_fm = (Po-fmj ^-fmi QI-fm = 0)> into the 
self-consistent mean-field equations, we find once again the fixed point is undetermined. 
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The equations for p c a _-pM an d s o-fm become dependent. The relation between p° a _-pM an d 
s o,-fm defines a line of fixed points joining [b + c + d, b, c, d] and </> q _fm [b + c + d,b, c, d] . 
The relation between p£_ FM and s^-fm is given by the condition 



(c + d)b 
led 



{Pa-FM ) ^ - 4 ^^<-FM + (1 - <-FM) 2 = . (68) 



Similar considerations hold for the other transitions occurring in the model. For 
instance, at the transition towards the c-AF phase, c c = a + b + d, one can identify a 
line of fixed points of the form </£af = {Pc-af = 0' S c-AF' 9c-Af)> where s^_ AF and <j£-af are 
constrained to fulfill a given condition similarly to (I68[) for the a-FM transition. 

For the disordered points lying on the surfaces a + d = c + b or a + c = d + b in the PM 
phase we have spm = 0. The free-energy at these points computed with the single vertex 
tree, eq. (|26p . is the same as the free-energy obtained using the plaquette model and 
coincides with the exact result of the 2D model [27]. This can be understood in terms 
of the transformation of eq. (139|) . which maps these particular points in the PM phase of 
the phase diagram into the frozen FM phase of the six- vertex model (see next section). 
Since the particular structurte of the graph is irrelevant for the FM (frozen) phase of the 
six-vertex model, the free-energy on the Bethe lattice coincides with the exact result on 
the square lattice in 2D. 



4.5 Duality in the eight- vertex model 

Both the tree of single vertices and the tree of plaquettes yield the same critical planes, 
given by the condition |A§| = 1, which also coincides with the exact result in 2D. This 
is surprising, since in general the location of the critical points depends on the character- 
istics of the lattice. This result can be interpreted in terms of a duality transformation 
connecting the disordered phase with a ~ b ~ c ~ d and the ordered phase a>i,c, d. 
Such transformation is independent on the particular structure of the lattice and only 
depends upon the connectivity. For the 2D model on the square lattice the duality is 
discussed in Sec. 10.2 of [27] , and it can be easily generalized to any graph of connectivity 
four. It connects the partition function of the model with fugacities a, b, c, d to the one 
with fugacities a',b',c' and d', under the mapping (I59"j) . The transformation (1591) is an 
involution and one may express a, b, c, d in terms of a', b', c', d' exactly in the same way. 

At the level of our computation with the plaquette one can note that T[a', b', c', d'\ = 
u[a,b,c,d] or, similarly, T[a,b,c,d] = v\d ',c',d']. Moreover, the duality holds for the 
free-energies /pmK, b', c', d'} =f a _FM[ a ,b,c,d]. Then, one can map one solution into the 
other, 

s PM [a,b, c, d] 



^ PM [a,6, c, d] -^™[a,b,c,d] 
^™[a,6, c, d] +^[a,b,c,d] 

■Sq-FM ~ Pa- F M[a',b',C?,d'} 

1 + p a -FM[a',b',c',d'} 

</a fm k, v, c ', d'\ - v, c ', d'] 



(69) 



^"+ M [a', V, C, d'} + i)X- M [a', V,d, d'] ' 
and vice versa, as well as the thermodynamic quantities. The transition point can be 
recognized as the fixed point of the transformation (139|) which is consistent with the exact 
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result in 2D, a c = b + c + d. Thanks to the mapping of eq. fl69l) the infinite temperature 
solution 0pm [a, a, a, a] can be mapped into the completely ordered state a -FM[a', 0, 0, 0]. 
However, at the transition point <ppu[b + c + d, b, c, d] ^ q -fm [b + c + d, b, c, d] and the line 
described by eq. (|68p connects the two solutions. The same duality holds for the solution 
of the single vertex Bethe lattice where the free-energy of the PM phase can be mapped 
into the free-energy of the completely frozen FM solution under the mapping given by 
eq. U32D. 

The fixed point of (|39|) gives the transition point for the a-FM transition. This is 
enough to recover the transition planes relative to the other ordered phases using the 
symmetries that connect the different types of ordered phases, as discussed in [27J. 

4.6 The sixteen-vertex model on a tree of vertices 

In the previous sections we compared our mean-field approaches with the exact results 
of the six- and eight- vertex models in 2D, showing that the BP computation allows one 
to describe all the expected phases of the models and yields remarkably accurate results 
compared to the 2D case. We now focus on the mean-field analysis of the complete 
sixteen-vertex which includes also charge ±2 defects, see Fig. HI 

4.6.1 Recursion relations and fixed points 

We call ip^ 3 the probability for the root vertex i on a rooted tree with missing edge 
(i a jP) of type ( with ( e xl 6 { v ij v 2, v ie}- Similarly to the analysis of the six- and 
eight- vertex models, we introduce the cavity probabilities ijjf = if)^ (+1) to find a 
positive spin on the missing edge of the rooted tree; these probabilities satisfy 

€ = €r r + + + + + <T jU + <7 r + <7 jU , 

# = €r jd + + €r jd + €r f + <: jd + + + <r jd , (70) 

4 = + + + + + + + <^ , 
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Along the same line of reasoning outlined in Sec. I4.3[ we derive the set of self-consistent 
equations for the probabilities defined above: 

i/j u = ij) u [a, b, c, d, e, ij; u , ifj d , ip r ] 

= —\a ^V> r + b (1 - i) l U v (l -ip r )+c(l- ip u )(l - i/> l U r + d ^(1 - ip u )(l - V) 
z u L 

+e - ip u )ip r + ip l ip u (l - + (1 - ip l )^ r + (1 - ^')(1 - - ip r )^j 

^ l = i> l [a,b,c,d,e,rA d A l A r ] 

= \\a ^VV" + 6(1- if> d W(l - ip u ) + c ip d (l - -i[) u ) + d(l- i) d )(l - ifj l )ij u 

+ e (Vy (i - + - *i> l )if> u + (i - ^ d )^v u + (i - - ^'x* - 

^ = y d [ a , b, c, d, e, V", V', V] 

= J_ L + ^> d (i - ^) + c (i - v r )(i - i> d W + d i> r {i - i) d ){\ - ip) 

z - 

+e ((1 - ip r )ip d ip l + ip r (l - ip d )if} 1 + %p r %p d (l - ip 1 ) + (1 - ^ r )(l - ^ d )(l - ^ 
^ = q u [ a ,b,c, d,e,^ u ,^ d ,^,*p r } 

aip u i/j r ?p d + 6(1 - */j u )ip r (l - */j d ) + ap u (l - i/j r )(l - i) d ) + d(l - i/j u )(1 - ip r )i/j d 
+e(V(l -ip r )*P d + *P U {^ -ip r )ip d + ip u {l -ip r )ip d + (1 -ip u ){l -ip r ){l -ip d )^j 

n (71) 

where z a are normalization constantsjj In order to describe AF phases, these equations 
must be studied on bipartite graphs as in eq. ff2T]) . 

For simplicity, here and in the following we consider the case in which the vertices 
Vg,...,Vi6 have the same weight e; arrow reversal symmetry holds for the statistical 
weights of all other vertices as well. However these assumptions could be easily released. 
The fixed points of eqs. (ITTj) on one of the two sub- lattices, say A\, are 



^PM z 


~ ^PM — 








V'a-FM 


= ^a-FM 


= CfM = V'a-FM = 


h + (a, b, c, d, e) , 






= ^6-FM 


= h + (b, a, c, d, e) , 


i'b-FM = ^Pb-FM 


= h_(b, a, c, d, e) 


C-AF 


= AF ~- 


= h + (c, b, a, d, e) , 


^c-AF = ^c-AF = 


= h-(c, b, a, d, e) 


V'd-AF 


= ^d-AF 


= h + (d, b, c, a, e) , 


V'd-AF = VVAF 


= h-(d, b, c, a, e) 



(72) 



with 



h±(w 1 ,w 2 ,w 3 ,w i ,w 5 ) 



± 



U>1 — iu 2 — If 3 — — 2w 5 
2\/ (wi -W2-W3- W4) 2 - Awl 



6 Thc z a differ from the ones defined in eq. (flT))) . due to the contributions from e 7^ 
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4.6.2 Free-energy, stability and order parameters 

In presence of vertices of kind e, the vertex partition function reads 

Z v [ip l ,ip r ,ip u ,ip d ] = a ifj l ^ u ^ r ip d + (1 -ip u )(l -ip r ){l -ip d ) 

+ 6(1- ip l )^ u (l - ip r )ip d + ^(1 - ip u )ip r (l - ip d ) 

+ c (1 - ip u ){l - ip l )ip r if} d + ^V(l - ^ r )(l - ^ d ) 

+ d - ip u ){l - ip r )ip d + (1 - ip l )ip u ip r {l - ip d ) 

+ e - tfj u )ilj r tp d + (1 - i/> l )if> u (l - - i> d ) 

+ - + (1 - ^)(! - ^ n )^ r (l - 

+ (1 - ^ l )i) u i) r i) d + ^(1 - ^ tt )(i - VOC 1 - 

+ (1 - ^)(1 - ^ n )(l - ^ r )^ d + - 
which in the limit e — >■ gives back eq. (123]) . The free-energy reads 

0f[a,b,c,d,e, i^i,^] = -^( ln ^['0i] + In ^h/> 2 ] + 

- In Z (ir> ^] - In Z (ir> VI] - In Z {ud) ^] - In Z {ud) fo$, < 



(73) 



(74) 



where Z^ r ) and Z( ud ) are defined in eqs. (124|) . The free-energies of the different phases is 
obtained evaluating eq. ( !74|) at the fixed points given by eq. (172]) : 

p/pm = P7 K o, c, d, e, VpmJ 



-In 

/3/a-FM = pYk b, c, d, e, Wfm] = - In 
Pfb-FM = 13 f [a, b, c, d, e, t/>&-fm] = - In 
/9/o-af = P7[a, 6, c, d, e, V^-af] = - In 
/3/d-AF = /9/[a, 6, c, d, e, Waf] = - ln 



2e 5 



a — 



c — 



d 



-a + 64 


-c + 


2e 2 




a — 6 + c 


+ dJ 


2e 2 




a — c + 6 


f d. 


2e 2 




a — d + b 


+ cJ 



(75) 



m a-FM 



The order parameters are defined as in Sec. I4.3.3l with the addition of the contributions 
from the new vertices. For instance, the FM order parameter m a _FM reads: 

— \a U l ip u ip r ip d - (1 - ip u ){l - - 4> r )(l - ij> d j) 
Z v l\ J 

+| ( - A - i>u - - i>i + 2^ r ^ n + 2^# r + 2ip u ipi + 2ip d if;i (76) 
+2^ d ip u + 2ip d ip r - 2^iijj r ijj u - 2ip d ip r ip u - 2ip d ipi^ r - 2VW#« 
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where the normalization is now the one defined in eq. fj73|) . The other order parameters 
characterizing the other ordered phases have a similar expression. The PM solution 
V'pm i s the same as for the eight-vertex model and yields vanishing magnetization for 
all the order parameters. Differently from the six- and the eight-vertex model, where 
the order parameters jump discontinuously to one at the phase transition on the tree of 
single vertices, when e 7^ any ordered FM or AF solution is associated to a continuous 
transition. For instance, for the a-FM phase transition, plugging i/Vfm from eq. (T72|) 
into eq. (f76i) one obtains: 



m a-FM 



^(a-b-c-d) 2 - 4e 2 
a — b — c — d 



(77) 



The expansion close to the critical point a c = b + c + d + 2e yields the expected mean- field 
classical exponent (3 = 1/2: 



m a-FM 



V a - a c 



+ 



( \ 3 / 2 

[a - a c ) 



(7£ 



In the limit e — ?• one recovers from eq. ( 177|) the value m a _ F M = 1, characteristic of the 
frozen (on the tree of single vertices) ordered phase in the eight-vertex model, due to a 
divergence of the coefficient of the singular term y/a — a c . The magnetizations tt^fm, 
m c _AF and rrid-AF are zero in the a-FM phase. Analogous results hold for the other phase 
transitions and the corresponding order parameters; these results can be obtained by 
simply exchanging the parameter a and the opportune vertex weight in eqs. fl77|) and 

(HE}. 

We study the stability properties of the PM phase by analyzing the eigenvalues of the 
stability matrix introduced in eq. (133]) for e > 0. The eigenvalues associated to the PM 
solution i/'pm are 



PM 



E 



PM 



3a — b — c — d 

a + b + c + d + Ae ' 

a + b — 3c + d 
a + b + c + d + Ae ' 



E. 



PM 



-a 



3b 



d 



E 



PM 



a + 6 + c + d + 4e 

a + b + c — 3d 
a + 6 + c + d + 4e ' 



(79) 



with the same eigenvectors as in Sec. 14.3.41 Overall the stability is controlled by the 
condition 



1 + E 3 PM )(1 + E 



PM> 
4 ) 



[1 -£f M )(l 



E™) 



with 



;i+£ 3 PM )(i 



A sv 

^16 



E™) 



E™) 



f (l-£f M )(l 
c 2 - d 2 + 2(a + b - c - d)e 



^16 1 



< 1 



2[cd + ab + e(a + b + c + d + 2e)] v ' 

This generalised anisotropy parameter implies that, in the presence of a non-vanishing 
vertex weight e, all the transition lines are shifted by 2e with respect to the value obtained 
for the eight-vertex model. This should be compared to the conjectured value for the 
2D model, A 16 , given in eq. f[T41) . and with the numerical results in 2D on the square 
lattice. The BP approach on the tree of vertices gives the correct qualitative behavior, 
meaning an increase of the PM phase with increasing e, but it does not coincide with the 
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numerical results. The conjectured Ai6 given in eq. (I14p is closer to the numerics (this is 
not surprising, since its functional form was guessed from the analysis of the numerical 
results). Unfortunately, there is no exact result in 2D for the sixteen- vertex model to 
compare with and we are not able to draw more precise conclusions on the anisotropy 
parameter 



4.7 The sixteen- vertex model on the tree of plaquettes 

The procedure described in the previous sections can be easily extended to analyze the 
sixteen-vertex model on the tree of plaquettes. We use the definition 



w si,s2,s3,s4 \ a i ^> c > e) 



a {I + S1S2S3S4) + 6'(sis 3 + s 2 s 4 ) + c'(sis 4 + s 2 s 3 ) 

(82] 

1 p 

+ d'(s 1 s 2 + S3S4) 



|(1 - S1S2S3S4) 



where a', b', d and d! given in eq. (155)) . Then, eqs. ( I4"3"j) and (JHJ yield the fixed-points, 
eq. (|45p gives the free-energies of the different phases and eqs. 047|) determine the order 
parameters of the sixteen-vertex model. As the expressions are rather long and involved 
we prefer not to write them down explicitly here. 



4.7.1 Paramagnetic phase 

The PM solution is of the type PM = </>p M = (f) l PM = PM = 0p M = (p PM = 
0, spm, qpM = 0) and spm is the solution of the equation: 

(a + b + c + d) 4 (x 2 - y 2 ) A 4 4m + ^3 4 M + ^2 sp M + Ai s PM + A =0 (83) 

with 

X {x,y,z,u) = (1 + uf + z 2 , 

(1 + u) A + z 4 - 2(1 + z 2 - u 2 )(x 2 + y 2 ) 



X!(x,y,z,u) 



x 2 — y 2 



X 2 (x, y, z, u) = -12u , (84) 

{l-uf + -2(1 + z 2 -u 2 ){x 2 + y 2 ) 

A 3 (x, y, z, u) = = -Ai(x, y, z, -u) , 

x l — y l 

Xi(x,y,z,u) = -[(1 -u) 2 + z 2 ] = -X (x,y,z,-u) . 

The variables x, y and z were already defined in eq. ( 1511) and u = 4e/(a + b + c + d). This 
equation reduces to that of the eight- vertex model, eq. (|5"3"1) . when u = 0. The presence of a 
non-zero value of u implies a much more complicated dependence on the parameters as the 
roots of eq. (183]) are more involved. We do not solve this equation explicitly; one can show 



7 Actually there is no guarantee of the existence of such a single Ai6 parameter for the sixteen-vertex 
model. 
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55) 



that, as for the eight vertex model, the limit of infinite temperature a = b = c = 
corresponds to the trivial solution spm = 0. At linear order in u one obtains 

1-Vr 8(x 2 - y 2 ) 

SPM ~ 1 + Vf + (1 + + Ax 2 + z 2 )(-l + Ay 2 + z 2 ) X 

(1 + 4(a: 4 - Qx 2 y 2 + y') - 4(x 2 + y 2 )z 2 - z') 

. . u + U[u ) , 

+ Ay 2 + z 2 + V-l + 4x 2 + z 2 ) 2 

where T is the function defined in eq. (154]) which characterizes the PM phase when u = 0. 
4.7.2 Ferromagnetic phase 

Solving numerically eqs. f l44"|) for the sixteen-vertex model one finds that the a-FM fixed 
point takes the general form p u = p l = p r = p d = p a_FM ; s u = s l = s r = s d = s a_FM , 
q u = q l = —q r = —q d = g«-FM w j^] 1 ^a-FM gjjg^iy different from zero unless c = d or 
e = 0. This implies that the corresponding eigenvector has a non-vanishing component 
also along the "direction" of q. However, g a_FM is very small and influences very little 
the transition point, while the instability is mainly driven by the FM component which 
gives a non zero value of p a ~ FM . Therefore, while in the numerical study we used the 
exact expression of the eigenvector, in the analytic study of some asymptotic behaviors 
we disregarded the components contributing to the non- zero value of g a_FM . Accordingly, 
in order to study the stability we consider the simplified condition: 



dffi 4> PM dp a 4>pm 

a=u,l,a,r a=u,l,d,r 

where /3 can be u, d, I or r (the result does not depend on 0). This condition translates 
into the equation: 

A 3 s 3 c + A 2 s 2 c + Ai s c + Ac = (87) 

with 

A (x, y,z,u) = -(1 + u) A + x(l + u) 3 + 2x(x 2 + y 2 + 2uy 2 ) + 

2(x(l + u) 2 + x 2 (l + u + x) + (1 + u- x)y 2 )z + (1 + u)(l + u + x)z 2 , 
Ai(x, y, z, u) = 4x 3 (l + z) + (y 2 - x 2 )(-l + Au + (1 + u) 2 + (z - l) 2 - 2z) + 

2x(2y 2 (-l + z) + z(-u 2 + (1 + z) 2 )) , 
\ 2 {x, y, z, u) = ((« - l) 2 - z 2 )z 2 + 2(x 2 + y 2 )(-(l - z) 2 + u 2 - z(l + u)) + x(l - u) 3 

+x (-1 + 2x 2 (l + z) - 2y 2 {-l + z + 2u) + (1 + zf - uz{4 + z- 2u)) 
~\ 3 (x,y,z,u)= (-x 2 + y 2 )((-l + u) 2 + z 2 ) . 

(88) 

Altogether, eqs. fl83|) and flHTj) give the value of s c and the critical point a^f, when the 
other four parameters are fixed. The 6-FM phase has the same properties as the a-FM one 
and the solution is of the form: p u = —p = —p r = p d = p b_FM ; s u = s l = s r = s d = s 6_FM , 
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q l = q d = —q u = —q r = g fe_FM with q b ~ FM very small but different from zero unless c = d 
or e = 0. If one exchanges the weights 6 with a while keeping fixed the other fugacities 
the quantities above are precisely the same as for the a-FM phase. 



4.7.3 Ferromagnetic transition in the limit a, 6 ^> c,d,e 

Let us focus on the FM transition at large values of a for the model on the tree of 
plaquettes. In order to extract some analytic behavior we will focus now on the limit in 
which one of the four other fugacities is large and, in particular, in the regime a, 6 3> c,d,e 
or a, c > c, d, e. This situation corresponds to the cases where two types of vertices have 
an energy much smaller than the others. In order to study the large b behavior, on the 
basis of the numerical results we take a to be of the form — b + r ab , where r a b(c, d, e) 
is some (small) function of the remaining parameters. In this limit one has x ~ 
y ~ z ~ 1 — ^ and u ~ y. In the infinite b limit eq. (183]) admits as the only 
acceptable solution spm = s c = 0, which does not depend on the parameters. This is also 
confirmed numerically, where one sees that in the limit a, b 3> c,d,e the value of s c at 
the transition point goes to zero. The same result is obtained in the eight-vertex model 
in the same limit. We then linearize eq. (j83p in s c and we expand the solution at leading 
order in 1/6. In conclusion, we find 

c " 8(c + d+2e)b ■ { ' 

Plugging this result into eq. (I8"7|) and taking the leading order in 1/6 one obtains r ab (c, d, e) - 
c + d + 2e, which implies 

a P l = a s c v = b + c + d + 2e for a, b > c, d, e, (90) 

as for the single vertex tree. In the left panel of Fig. [TBI we show the comparison between 
the exact a-FM transition point for the plaquette (solid lines) and the one predicted by 
the single vertex (dashed lines), for different values of the vertex weights c,d and e. In 
the limit a,b ^> c,d,e the two results turn out to be remarkably close. 



4.7.4 Ferromagnetic transition in the limit a,c^> b,d,e 

This limit is more involved and shows a different behavior from the one found above for 
both a and b large. Supported by the numerical results we assume that the critical a c 
scales as a% 1 = c + r ac , with r ac (b, d, e) some function of the other parameters. In this 
limit x ~ \ — 3&+ 4~ r ° c , V — \— b+3 ^+ Tae , z ~ b ~ d + T °>c and u ~ y. Plugging these results 
in eq. (183)) one sees that in limit of large c, differently from above, s c does not go to zero 
but it converges towards a finite value which has to satisfy the relation 

b - d - r ac + 4(6 + d + 4e)s c - 4(6 + d - Ae)s 3 c + (-b + d + r ac )s 4 c = . (91) 

Even though u goes to zero with c, the solution to this equation is not the same as the 
one of the eight-vertex model in the limit of large c and = c + r ac , since an explicit 
dependence on e remains. The large c limit in eq. (|87|) yields 

(-3b-5d-8e + 3T ac ) + (b-d + 3T ac )s c + {b-d-8e + 3T ac )s 2 c + {b-d-T ac )s 3 c = . (92) 
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Figure 18: Left panel: critical value a?c versus b on the tree of plaquettes for different values 
of the other parameters c, d and e (solid lines) compared to a^ 1 ' = b + c + d + 2e found on the 
single vertex tree (dashed lines). The plot shows that in the large b limit the FM transition 
is well described by af. Right panel: critical value a£ versus c on the tree of plaquette for 
different b, d, and e (solid lines) compared to the asymptotic behavior following from eqs. (j93[) 
and (|94p (dashed lines). The variation of the prefactor multiplying e shows that a linear fit in 
all fugacities (as for the single vertex approximation) does not work in this regime. Note that 
such non linear effects can also be seen in the left panel, in the small b regime, where c is the 
dominating weight. The agreement with the asymptotic limit improves for increasing c. This 
proves that the behavior a£ = c + T ac (b, d, e) is justified. 



We can now use eq. f[9~Tj) to determine r ac , obtaining 

16es c 4(b + d)s c , . 

Tac = b-d+- 1+ \,i C 93 

i -s 2 c i + s 2 c 

which substituted in fl92|) gives an equation for s c : 

2b + d-5e — — + (2b + d - 6e)s c + (-d + e)s 2 c — ; + / = . (94) 

— l + s c 1 + s c 1 + si 

Both equations (1931) and (1941) are linear in the parameters b, d and e; instead, the function 
r ac will generically show a non linear behavior, brought about by the non trivial value of 
s c (b, d, e), which is the asymptotic result in the large c limit. In the right panel of Fig. [TH] 
we show the critical lines for the a-FM transition obtained for the plaquette (solid lines) 
and from the asymptotic expressions extracted from Eqs. fl93l and f[94"j) for different values 
of b, d and e (dashed lines). 

Let us finally mention that another result can be obtained when b = c = d = and 
e > 0. In this limit we take a v ^ = r ae e and solving eqs. fl8"3"|) and fl8T|) one finds r ae ~ 2.34 
and s c ~ 0.144. Apart from this "extreme" case we did not explore the large e regime. 

4.7.5 Antiferromagnetic phases and transitions 

The AF phases are characterized by similar solutions. In particular the phase dominated 
by vertices of type c is described by the fixed point p u = p l = —p r = —p d = p c ~ AF , 
s u = s l = s r = s d = s c ~ AF , and q u = q l = q r = q d = q c ~ AF , with p c ~ AF very small but 
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different from zero, unless a = b or e = 0. As we discussed for the FM transition, in 
order to identify the point of instability one can disregard the component along p as it is 
negligible for the AF transitions and study the quantity 

= E Sri = E =1. 

^— ' d<7>? <4 PM ^-^ dq a <7!)p M 

a=u,l,d,r r3 WM a=u,l,d,r H VPM 

where /3 can be it, d, I or r. Equation (1951) can be obtained from eq. (IHTj) by exchanging 
a with c and 6 with d and sending s to — s. 

The same results discussed for the critical FM surfaces hold when one considers an AF 
transition (for instance the one at large c). In this case the critical value dP^ obtained with 
the plaquette coincides with the one of single vertex in the limit in which c,d 3> a,b,e 
while the case c, a ^> b,d,e shows more subtle non-linearities. This is in fact the case for 
the AF transition lines reported in Fig. [T9j 

4.8 The sixteen vertex model: phase diagram and discussion 

In Fig. [TjJ]we plot the projection of the transition lines and the resulting phase diagram 
of the sixteen- vertex model on the plane (a/c,b/c), for two fixed values of d/c = e/c, 
showing the a-FM, the b-FM, the c-AF and the PM phases. The transition lines within 
the plaquette approximation can be either extracted from the analysis of the stability 
matrix (149)) or one can directly detect the point where the self-consistent equations develop 
a symmetry breaking solution. The results obtained with the tree of plaquettes are shown 
with orange (d/c = e/c = 0.1) and red (d/c = e/c = 0.2) dotted lines; the transition lines 
obtained according to the guessed value of the anysotropy parameter given in eq. (TMj) 
are drawn with green (d/c = e/c = 0.1) and violet (d/c = e/c = 0.2) straight lines while 
the exact phase diagram of the six-vertex model for d = e = is shown with blue solid 
lines. The same code of color - green, violet and blue - is used to represent with squares, 
stars and triangles the results obtained by MC simulations on the 2D square lattice. In 
addition we also indicate with a black dot an extra point obtained by MC simulations 
for the c-AF transition with d/c = e/c = 0.05. The results obtained with the single 
vertex are not shown in Fig. [19] however the a-FM and b-FM transitions lines found with 
the plaquette are in very good agreement with those of the single vertex in the regime 
a, b ^> c, d, e. 

The phase diagram of the model defined on the tree made of single vertices and of 
plaquettes is in qualitative agreement with the MC simulations on the 2D square lattice. 
All the transitions are continuous when all possible vertices are present. 

The results on the tree of vertices predict a uniform shift of the critical lines, with 
respect to the eight-vertex model, by 2ej^| This implies, for instance, that the critical 
value of the a-FM transition occurs at a s c v = b + c + d + 2e. This does not reproduce 
exactly the numerical results for the 2D model on the square lattice. While for small 
values of d and e, one numerically finds a 2 c D ~ b + c + d + 3e (the results reported in |29j), 
for large values of these fugacities we see a deviation from the linear behavior. We do not 
have an analytic expression for the transition surfaces for the tree of plaquettes but they 
can be found from the numerical solution of the self-consistent equations with arbitrary 

8 This is obtained by comparing cq. ([79)) with the criticality condition \E a \ = 1. 
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Figure 19: Phase diagram of the sixteen-vertex model. The figure shows the projection of the 
transition surfaces on the plane of parameters (a/c, b/c), for two fixed values of d/c = e/c = 0.1 
and d/c = e/c = 0.2. Orange and red dotted lines represent the results obtained with the cavity 
method for the tree of plaquettes. The critical lines obtained with the tree of single vertices are 
not shown but they are in very good agreement with the results obtained with the plaquette 
for the ct-FM and 6-FM transitions in the regime a,b 3> c, d, e. Green and violet plain lines 
shows the proposed behavior of the transition lines predicted by eq. (|14p for the 2D model. 
Blue solid lines correspond to the exact phase diagram in the limiting case d = e = (six- vertex 
model). The same code of colors is used to indicate the transition points obtained with Monte 
Carlo simulations for the model defined on the square lattice. In particular green squares are 
for d/c = e/c = 0.1, the violet star for d/c = e/c = 0.2 and blue triangles for d = e = 0. We 
also indicate with a black dot a c-AF transition point obtained with d/c = e/c = 0.05. 

precisions and their asymptotic expressions can be obtained in some particular limits. In 
particular, when a,b 3> c,d,e for the FM transitions and when c,d 3> a,b,e for the AF 
ones, one recovers the results of the single vertex. It would be interesting to see whether 
asymptotically in the same regime a similar result holds for the 2D lattice. Other cases 
are generically accompanied by (small) non-linearities in the fugacities. Comparing our 
analytical approaches with the MC data, we conclude that for the FM transition one finds 
a s c v < < a 2 c D (where the superscript sv, pi and 2D stand for single vertex, plaquette 
and 2D model respectively). Similar results hold for the other transitions. 

4.9 Fluctuations in the symmetry- broken phases 

In this subsection we discuss the nature of thermal fluctuations in the symmetry broken 
phases of the sixteen-vertex model and we compare them to the ones found in the PM 
phase. For simplicity we focus on the a-FM phase only, since the observations reported in 
the following can be extended to the other ordered phases, as we show in the figures. As 
we already mentioned, the solution 4>2-fm within the plaquette approximation is charac- 
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terized by FM order signaled by a finite value of p a -FM and by a small but non- vanishing 
value of the AF field g a _FM- This field has not the same sign along all the directions, 
and overall it does not lead to a staggered magnetization. Moreover, it decreases upon 
increasing the FM order. Interestingly enough, such small AF field plays a crucial role 
for thermal fluctuations in the symmetry broken phase. 

In Fig. |20]we show some equilibrium MC configurations in the a-FM phase on the 2D 
lattice, for different values of the ratio c/d. The comparison between the three snapshots 




Figure 20: MC a-FM equilibrium configurations of the 2D sixteen-vertex model for different 
values of the ratio c/d, showing that fluctuations are anisotropic for c/d ^ 1. The color code 
is the one used in Figs. and [H blue is for vertices b, orange for c, green for d and red for e. 
FM vertices of type a which are colored in black (v\) and white (1*2) • The size of the system 
is L 2 = 50 2 . Left panel: a = 2.5, b = 0.3, c = 1.5 > d = e = 0.1. Linear fluctuations made 
of c vertices extend along the down-left ^ up-right diagonal direction. Central panel: a = 2.5, 
b = c = d = 0.5, e = 0.1. Fluctuations are homogeneous with no preferential direction. Right 
panel: a = 2.5, b = 0.3, d = 1.5 > c = e = 0.1. Linear fluctuations extend along the down-right 
^ up- left diagonal and are made of d vertices. 



clearly shows that vertices c and d compete to form fluctuations along opposite diagonal 
directions (in the left panel c > d while in the right panel c < d) and that isotropy 
is recovered when their statistical weights are the same (central panel). One can also 
find fluctuations of e vertices alone, which are generically unstructured, or of b vertices, 
that consist in straight vertical and horizontal lines of defects with the same probability. 
In the following we concentrate on diagonal defects which are the origin of anisotropic 
fluctuations, and which are related to the small AF field found in the mean-field solution 
of the model on the tree of plaquettes. 

In the a-FM phase, when c> d, fluctuations constituted of diagonal strings of defects 
of c vertices grow along the diagonal running from the up-right corner to the down-left 
one. They can be mainly of two types: either one-dimensional "stairs" of arrows going in 
the opposite direction with respect to the ordered background and extending between two 
vertices of type e; or thicker defects containing domains of FM vertices with the opposite 
magnetization with c vertices on the boundary (see Figs. [2D] and I2TT) . As can be seen 
from Fig. ETJ in both cases vertices of type v 5 are above vertices of type v q with respect 
to the diagonal axis that joins the two extremities of the elongated defect. Reversing all 
the spins of the background (this corresponds to the situation in which the opposite FM 
order is the dominating minimum) amounts to invert the role of and Vq. This implies 
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Figure 21: (Color online.) Anisotropic fluctuations in the a-FM phase of the 2D sixteen-vertex 
model with c > d. We color the vertices according to the code used in Figs. and HI We 
distinguish the two vertices connected by arrow reversal symmetry by drawing a white dot in 
the middle. Arrows carry the color of the vertices in the background if they are compatible with 
the ordered phase, or the color of the neighboring vertices if they have an opposite orientation. 
The first and second panel show two staircase lines of defects made by an alternating sequence 
of c vertices (orange dots) between two vertices of type e (red dots). In terms of arrows such 
defects can be viewed as a one dimensional domain. The staircases can also be seen as double 
diagonal string of vertices along the down-left ^ up-right direction with ^5 above vq in both 
cases. The third panel shows a thicker defect string. The AF vertices of type c with opposite 
orientation are now split apart and constitute the two boundaries of a domain containing FM 
vertices with opposite orientation with respect to the background. The two extremities of the 
domain are here sealed by two opposite ci-AF vertices (green dots). In the sixteen-vertex model 
more complicated structures can also be found. Note that if the symmetry were broken in favor 
of vertices vi one would find the same type of fluctuations, with v 5 and vq exchanged. 

that the breaking of reversal symmetry, accompanied by the selection of a given order, 
is also accompanied by its own characteristic fluctuations, and by the selection of one 
particular pattern of AF vertices among the two equivalent ones by arrow reversal. This 
does not imply an overall staggered magnetization since the defects strings correspond to 
patches of the two opposite AF configurations with equal probability (as in the first two 
upper panels in Fig. |2~T1) . When d > c one has the very same scenario, with fluctuations 
organized in the perpendicular direction and made of sequences of v-j and v$ vertices. A 
similar phenomenon leading to striped domains takes place also out-of-equilibrium for 
the coarsening dynamics in the a-FM phase |29j . 

Such anisotropic fluctuations cause the loss of the homogeneity of the solution 4>a-FM 
in the two directions. Isotropy is recovered in the case c = d, when there is no preferred 
axis, as shown by the MC snapshots (correspondingly, the field g a -FM vanishes in this 
case) . 

The AF field observed in the mean-field solution of the tree of plaquettes in the 
a-FM phase, is the signal of such fluctuations that are induced by the anisotropy and 
the symmetry breaking characteristic of the ordered phase. The configurations that 
contribute the most are those in which one vertex of type c is sitting in the down-right 
corner of the plaquette or the opposite vertex c is in the upper-left corner, together with 
other three vertices of type a (see Fig. I2"3"j) . 

At the level of the tree of single vertices such anisotropy can be detected by looking 
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Figure 22: Anisotropic fluctuations in the c-AF phase of the 2D sixteen- vertex model. The 
diagonal fluctuations are now driven by a vertices (for a > b). The first and the second panel 
show staircases made of v\ or V2 vertices and ending on e vertices. On the diagonal joining the 
two vertices e there can be both of them, while away from it, v\ and Vi occupy well defined 
positions which depend on the ordered phase of the background (equivalently, in the AF case 
the positions of v\ or V2 vertices depend on which of the two sub-lattices the diagonal strings 
belong). The third panel shows an example of a more spread domain, analogue to the one 
presented in the third panel in Fig. I2T1 with FM vertices on the boundary and AF ones in the 
interior. 




i>+- ip++ i>-+ i>++ 



Figure 23: (Color online.) Proposal for the anisotropic configurations of defects on the plaque- 
tte that contribute the most to create an effective (small) AF field q a -FM m the a-FM solution. 
Black dots represent FM vertices associated to the symmetry broken phase, orange and green 
dots are AF vertices, of type c and d, that correspond to fluctuations. 

at the derivatives dip 01 /dip^\^ a , which take different values depending on a and f3 
and which ultimately determine the susceptibility along certain paths. In particular, the 
"diagonal" derivatives d^" / dtp r and d^/dip 1 become equal only when c = d, while for c > 
d one has di/j u /dip r > dip u /dip 1 . Such derivatives are generically different along different 
directions also when evaluated in the PM solution signaling that also the disordered phase 
is anisotropic (apart from some particular choices of the parameters), as shown in Fig. [2H 
In this case, if a > b, c, d, e it will be more likely to have vertices of type a, albeit 
there is no spin reversal symmetry and a finite global magnetization since vertices V\ and 
V2 appear with the same frequency. Therefore, the diagonal fluctuations at the boundary 
and within these domains can be either the ones typical of domains of v± vertices or those 
found in domains of vi vertices. All AF patterns (and the AF boundaries in Fig. 123]) are 
equally probable and the position of v§ or t> 6 can be interchanged without any energy 
cost. 
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Figure 24: (Color online.) MC equilibrium PM configurations in the 2D sixteen vertex model 
with a = 2, b = 0.7, c = 1.2, d = e = 0.2 and L = 50. The left panel shows the lattice colored 
according to the vertex configurations and following the convention used in Figs. [3j 01 and [201 
where we colored v\ vertices in black and vi in white (note that these vertices are the ones with 
the largest weight). The right panel shows the same configuration where we colored the arrow 
sites, in blue if Sij = +1 and in red otherwise. The two figures demonstrate that the PM phase 
is anisotropic as well. For a > c > d it is more likely to have boundaries between the domains 
of positive and negative arrows along the diagonal that joins the upper-right corner with the 
down- left one. Such boundaries are formed by c vertices (in orange in the left panel). 

5 Summary and conclusions 

Let us start by briefly reviewing the main results found in this paper. 

Table [3] summarizes the phase transitions found in the six- and eight- vertex models on 
the tree of single vertices and plaquettes. The entries in blue correspond to the transitions 
which are of the same order on the Bethe lattice and on the 2D square lattice, while red 
entries correspond to the transitions that are of different type. 





6V single vertices 


8V single vertices 


6V plaquettes 


8V plaquettes 


PM-FMs 


Frozen-to-sPM 


Frozen-to-PM 


Frozen-to-sPM 


Discont-to-PM 


PM-AFs 


Frozen-to-sPM 


Frozen-to-PM 


Cont.-to-sPM 


Discont.-to-PM 



Table 3: Transitions in the six- and eight-vertex models on tree-like graphs. We call sPM the 
(critical) PM phases with a maximal eigenvalue of the stability matrix identical to one, which 
are reminiscent of the critical (Coulomb) PM phases in 2D. The comparison with the behavior 
of 2D models on the square lattice is encoded by the following color rule: the entries of the 
table are colored in blue whenever the transitions of the 2D model and on the Bethe lattice are 
of the same kind; contrarily, we use red when the transition of the models on the tree are not 
of the same type as the 2D one. We recall that in 2D the PM-AFs transition in the six-vertex 
model is of Kosterlitz-Thouless type while it is second order in the eight-vertex model. We 
stress that the plaquette model provides an important improvement compared to the single 
vertex approach as, for instance, it is able to capture the anisotropic fluctuations in the ordered 
phases. 

In the 2D six- vertex model the disordered phase is critical (a SL), the FMs are frozen 
and the AF phase has fluctuations. In the tree of single vertices the PM phase (that we 
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call sPM in the table) is similar to a critical phase in the sense that one of the eigenvalues 
of the stability matrix is identical to one, yielding a diverging susceptibilityH The FMs 
and AF phases are frozen. In the mean-field treatment with the plaquettes the pseudo- 
critical nature of the PM phase is maintained, the FM phases remain frozen but the AF 
phase is not, as in the 2D model. 

In the 2D eight-vertex model the disordered phase is no longer critical and the FM 
and AF phases are no longer frozen. Moreover the transitions are all continuous. In the 
trees of single vertices and plaquettes the PM phase is also a conventional one, meaning 
that the stability analysis does not lead to a critical eigenmode. As for the FM and AF 
phases, they are frozen on the tree of single vertices but they are not frozen on the tree 
of plaquettes. With the latter improved treatment we find discontinuous transitions that 
are also associated to instability, as in other KDP-like transitions. However, differently 
from the cases encountered before, the transitions are now towards non-frozen ordered 
phases. 

In conclusion, in these two (integrable) cases the plaquette version of the mean-field 
approach allows us to introduce fluctuations in the ordered phases, making the transitions 
softer and closer to the ones in the 2D case. 

The location of the transition lines for the six- and the eight-vertex models is in- 
dependent of the particular structure of the lattice (provided that it has coordination 
equal to four) and it is identified by the condition |A 8 | = 1. This shows that the critical 
planes are insensitive to the local geometry of the lattice and to the length of the loops, 
since they only depend on the local connectivity of the graph. This result can be un- 
derstood in terms of the duality transformation discussed in [27] which can be extended 
to the model on the tree. On the contrary, for the sixteen-vertex model we found that 
the location of critical points depends on the graph and a duality transformation as the 
one that holds for the eight-vertex model is not easily generalizable to the "complete" 
model for all values of the parameters. While in 2D the criticality condition |Ag| = 1 is 
defined by the singularity of the free energy, in our mean-field approach we recover it as 
a special combination of the eigenvalues of the stability matrix that captures the insta- 
bility towards the four possible ordered phases. We used the same procedure to extract 
a mean-field A^g that characterizes the critical surfaces of the sixteen-vertex model on 
the tree. Such A^g might provide a good estimate of the transition parameters of the 
spin-ice model on the pyrochlore lattice, the properties of which are well- described by a 
mean-field treatment |53j . 

In the finite dimensional problem the SL-FM transition becomes second order as soon 
as d > or e > 0. We find that the critical exponents depend on the particular values 
of the fugacities, although we are not able to determine their explicit dependence, due to 
the difficulty of the numerical analysis. A real-space RG method could be helpful to deal 
with this question precisely. The fact that all transitions in the sixteen-vertex model are 
continuous is confirmed by the calculations on the tree of single vertices and on the tree 
of plaquettes (cavity method). 

The results obtained with the cavity method on the Bethe lattice of single vertices 
generalize some of the calculations of Refs. [121 EH E2J [651 ESj for the pyroclore system. In 

9 Models defined on a tree cannot have a power-law decay of the correlation functions as in the critical 
phases of finite dimensional systems, due to the fact that the number of neighbors of a given site at a 
given distance grows exponentially with the distance itself. 
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particular, such treatment was used to study the so-called Kasteleyn transition in presence 
of a magnetic field [53], and to discuss the properties of multicritical points and KDP 
transitions found on the pyrochlore lattice [66]. Moreover in [65], disregarding vertices of 
type d, a critical transition temperature towards the a-FM phase was determined, which 
is in agreement with our a s c v if one imposes the same approximation. Notice that while 
in [531 ESI [66] the original motivation was the study of the pyroclore lattice, ultimately the 
mean-field theory that is developed is equivalent to the one that we used for the 2D lattice 
(when we considered the single vertex problem and with no distinction among the different 
direction thus not giving access to staggered ordered phases). This is because one ends 
up with a tree structure of vertices with connectivity four and spins (arrows) shared by 
two neighboring vertices. Contrary to these papers, our approach keeps track of the four 
different directions and allows us to study simultaneously all possible phases, including 
the AF ones that are relevant for artificial spin-ice samples in 2D [67]. More details about 
such experimental applications will be discussed elsewhere [31] . In addition, our approach 
allows us to access the non-homogeneity brought about by certain preferential paths. 
Indeed, as we discussed in Sec. 14.9} we found that both the PM and the ordered phases are 
anisotropic (except for some special combinations of the fugacities) and domains of arrows 
or vertices and thermal fluctuations of string defects develop along some preferential axes. 
The same anisotropy is the one that drives the coarsening dynamics out-of-equilibrium, 
as found in [2S], and should be easy to visualize experimentally. 
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A The CTMC algorithm 



In this Appendix we give some details on the implementation of the Continuous Time 
Monte Carlo algorithm that we used to study the phase diagram of the sixteen vertex 
model. 

We chose to use single spin updates. We can easily predict the time needed to flip 
an arrow. Suppose that the system is in state /i(to) at time to- The exact probability of 
leaving the state \i after At trials is 

(W(» /i)) A< (1 - W{pi fi)) = (WQjl /i)) Ai - (WQjl -> /i)) Am . 

In order to get an estimate for this quantity we have to generate a random number £ 
uniformly distributed between and 1. The latter corresponds to At trials if < £ < 
(^V) A * - (W^) At+1 then At + 1 < j^f- < At. It follows that the number of steps 
needed to flip an arrow should be computed by 

At = Int I —. ^ r I + 1 (A.r 



,ln(l-E?i^(/i^^ 




Equation flA.lj) shows that we need to know the transition probabilities for all possible 
arrow-flips at each step. This is the main difficulty for the implementation of the CTMC 
algorithm. A simple choice suggested by the detailed balance is 

W(jjl -> //) = g(fx -> //) (A.2) 

where we have split the transition probabilities into an edge-selection probability 

qQjl fM 1 ) = 1/2L 2 , V/ (A.3) 

and a flip-acceptance probability 

AU ^a I ) = { ***>{- PW) ~ EQi))} if £(//) ~ > 
[ 1 otherwise 

The transition probabilities defining the dynamics of the system can now be written as 

W{n -)• n 1 ) = ^min (l, e -^')~E(^ ( A 5 ) 

satisfying detailed balance and ergodicity. This transition probabilities are the same as 
for the FSMC algorithm, but we should notice that in this case the aim is to compute 
the time we have to wait before we do a flip instead of sampling the configurations 
of the system. This transition probabilities are expressed as a function of the state 
of the chosen spin /. We have to compute the probability to stay in the same state 
W (//—>■ (J,) — 1 — X)/=i —> fi 1 ) to obtain At using eq. flA.lj) and then we need to 
know every possible energy change a single flip can produce. In the sixteen-vertex model 
there is a finite number of such possible processes (and then a finite number of possible 
transition probabilities) independently of the system size. This procedure can then be 
applied by making a list of all the arrows classified by their state, noted from now on I 
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and denned by its neighborhood (i.e. the type of its two adjacent vertices). Since each 
vertex can take sixteen different configurations, there are 8x8 such states for vertical and 
horizontal arrows, so a total of 64 states for each type of arrow. Following the original 
name of this method \5U\ this algorithm is a 256-fold way. The transition probability 
of the process \i — > /i 1 only depends on the state I of the J-th arrow before the flip. This 
can be clearly seen by rewriting 



exp [-f3(E^D - En)] = exp 



-He 



+ E 



vtP -E[V£]-E[V£l 



here E 



AD 



is the energy of the first adjacent vertex of the 7-th arrow after the flip 

from the state \i. To know the type of the neighboring vertices V\j and V 2 j at state /i is 
equivalent to know the state of the concerned arrow before the flip (the vertex types of 
the neighboring vertices after the flip are determined by the vertex types before the flip): 
the energy change after a flip depends only in its initial state. We define 



1 

2iV 



mm 



(l,e-^) 



where E\ is the energy difference after flipping an arrow in state /. It is useful for the 
implementation to note that we can compute At by counting the number of arrows 
occupying each one of the different possible states at each step. We substitute the latter 
equation by 



Q 



2N 

E 

1=1 



W(n -> \i 



(/)> 



256 

E 

1=1 



9i Pi 



(A.6) 



where g\ is the number of arrows in state /. We then need to keep record of the state of 
every arrow on a list at each step. After a transition this list must be updated. 



59 



References 

[1] S. T. Bramwell et al., Spin Ice, in Frustrated spin systems, World Scientific, 
2004. 

L. Balents, Nature 464, 199 (2010). 

M. J. Harris, S. T. Bramwell, D. F. Mcmorrow, T. Zeiske, and K. W. 
Godfrey, Phys. Rev. Lett. 79, 2554 (1997). 

S. T. Bramwell and M. J. Gingras, Science (New York, N.Y.) 294, 1495 
(2001). 

M. J. P. Gingras, Spin ice, Springer- Verlag, 2010. 

B. C. den Hertog and M. J. P. Gingras, Phys. Rev. Lett. 84, 3430 (2000). 
S. V. Isakov, R. Moessner, and S. Sondhi, Phys. Rev. Lett. 95, 217201 (2005). 
J. D. Bernal and R. H. Fowler, J. Chem. Phys. 1, 515 (1933). 
L. Pauling, J. Chem. Phys. 57, 2680 (1935). 

W. F. GlAUQUE and J. W. Stout, J. Am. Chem. Soc. 58, 1144 (1936). 

A. P. Ramirez, A. Hayashi, R. J. Cava, R. S. Siddharthan, and B. S. 
Shastry, Nature 399, 333 (1999). 

J. C. Slater, J. Chem. Phys. 9, 16 (1941). 

S. V. Isakov, K. Gregor, R. Moessner, and S. Sondhi, Phys. Rev. Lett. 93, 
167204 (2004). 

R. W. Youngblood, J. D. Axe, and B. M. McCoy, Phys. Rev. B 21, 5212 
(1980). 

R. W. Youngblood and J. D. Axe, Phys. Rev. B 23, 232 (1981). 

T. Fennel, P. P. Deen, A. R. Wildes, Schmalz, D. Prabhakaran, A. T. 
Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell, Science 
326, 415 (2009). 

F. H. Stillinger and M. A. Cotter, J. Chem. Phys. 58, 2532 (1973). 

D. A. Huse, W. Krauth, R. Moessener, and S. L. Sondhi, Phys. Rev. Lett. 
91, 167004 (2004). 

C. L. Henley, Ann. Rev. Cond. Matt. 1 (2010). 

C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 3 (2008). 



60 



[21] D. Morris, D. Tennant, S. Grigera, B. Klemke, C. Castelnovo, 

R. MOESSNER, C. CZTERNASTY, M. MEISSNER, K. RULE, J. HOFFMANN, 

K. Kiefer, S. Gerischer, D. Slobinsky, and R. Perry, Science 326, 411 
(2009). 

[22] R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, 
M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, 
Nature 439, 303 (2006). 

[23] E. H. Lieb, Phys. Rev. Lett. 18, 692 (1967). 

[24] E. H. Lieb, Phys. Rev. Lett. 18, 1046 (1967). 

[25] E. H. Lieb, Phys. Rev. Lett. 19, 108 (1967). 

[26] B. Sutherland, Phys. Rev. Lett. 19, 103 (1967). 

[27] R. J. Baxter, Exactly Solved Models in Statistical Mechanics, Dover, 1982. 

[28] E. H. Lieb and F. Y. Wu, Two dimensional ferroelectric models, in Phase Tran- 
sitions and Critical Phenomena, edited by C. Domb and M. Green, Academic 
Press, 1972. 

[29] D. Levis and L. F. Cugliandolo, EPL 97, 30002 (2012). 

[30] G. M. Wysin, W. A. Moura-Melo, L. A. S. Mol, and A. R. Pereira, 
larXiv:1208.6557l 2012. 

[31] D. Levis, L. F. Cugliandolo, L. Foini, and M. Tarzia, (2012). 

[32] F. Rys, Helv. Phys. Acta. 36, 537 (1963). 

[33] C. Nisoli, J. Li, X. Ke, D. Garand, P. Schiffer, and V. Crespi, Phys. Rev. 
Lett. 105, 1 (2010). 

[34] J. P. Morgan and C. H. Marrows, Private Communication (2012). 

[35] C. M. Yung and M. T. Batchelor, Nucl. Phys. B 435, 430 (1995). 

[36] V. Korepin, Soviet Physics Doklady 27, 612 (1982). 

[37] V. Korepin and P. Zinn- Justin, J. Phys. A 33, 7053 (2000). 

[38] R. J. Baxter, Phys. Rev. Lett. 26, 832 (1971). 

[39] R. J. Baxter, Ann. Phys. 70, 193 (1972). 

[40] M. Suzuki, Prog. Theor. Phys. 51, 1992 (1974). 

[41] B. Sutherland, J. Math. Phys. 11, 3183 (1970). 

[42] C. Fan and F. Y. Wu, Phys. Rev. B 2, 723 (1970). 



61 



[43] P. W. Kasteleyn, J. Math. Phys. 4, 287 (1963). 

[44] F. Y. Wu, Phys. Rev. Lett. 22, 1174 (1969). 

[45] M. P. Bellon, J. M. Maillard, and C. M. Viallet, Phys. Lett. B 281, 315 
(1992). 

[46] F. Y. Wu, Phys. Rev. Lett. 24, 1476 (1970). 

[47] G. T. Barkema and M. E. J. Newman, Phys. Rev. E 57, 1155 (1998). 

[48] O. Syljuasen and M. Zvonarev, Phys. Rev. E 70, 1 (2004). 

[49] H. G. Evertz, G. Lana, and M. Marcu, Phys. Rev. Lett. 70, 875 (1993). 

[50] D. Allison and N. Reshetikhin, Annales de Vlnstitut Fourier 55 (2005). 

[51] M. Suzuki, Phys. Rev. B 31, 2957 (1985). 

[52] M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976). 

[53] L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, 
Phys. Rev. Lett. 100, 067207 (2008). 

[54] R. G. Melko and M. J. P. Gingras, J. Phys. A 16, R1277 (2004). 

[55] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comp. Phys. 17, 10 (1975). 

[56] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical 
Mechanics, Oxford Univ. Press, San Diego, 1999. 

[57] H. K. Janssen, B. Schaub, and B. Schmittmann, Z. Phys. B 73, 539 (1989). 

[58] A. Jaster, J. Mainville, L. Schulke, and B. Zheng, J. Phys. A 32, 1395 
(1999). 

[59] E. V. Albano, M. A. Bab, G. Baglietto, R. A. Borzi, T. S. Grigera, E. S. 
Loscar, D. E. Rodriguez, M. L. R. Puzzo, and G. P. Saracco, Rep. Prog. 
Phys. 74, 026501 (2011). 

[60] H. A. Bethe, Proc. Roy. Soc. London A 150, 552 (1935). 

[61] M. Mezard and A. Montanari, Information, Physics and Computation, Oxford 
University Press, 2009. 

[62] S. Yoshida, K. Nemoto, and K. Wada, J. Phys. Soc. Japan 73, 1619 (2004). 

[63] L. Zdeborova and M. Mezard, Phys. Rev. Lett. 101, 078702 (2008). 

[64] L. Benguigui, Phys. Rev. B 16, 1266 (1977). 

[65] L. D. C. Jaubert, Topological constraints and defects in spin-ice, PhD thesis, 
2008. 



62 



L. Jaubert, J. Chalker, P. Holdsworth, and R. Moessner, Phys. Rev. Lett. 
105, 87201 (2010). 

D. Levis, Two dimensional spin-ice and the sixteen-vertex model, PhD thesis, 2012. 



63 



